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1. _ Introduction 

Computer  molecular  dynamics  (CMD)  is  a  powerful  technique  for 
calculating  equilibrium  and  dynamical  properties  of  solids.  The 
method  consists  of  calculating  the  classical  trajectories  of  several 
hundreds  of  particles  interacting  through  a  known  potential  function, 
by  numerical  integration  of  the  equations  of  motion. 

The  equations  of  motion  which  govern  the  time  evolution  of  the 
nuclear  coordinates  are  a  set  of  coupled,  non-linear  ordinary  dif¬ 
ferential  equations.  The  time  evolution  of  the  phase  coordinates  is 
called  the  trajectory,  and  the  CMD  technique  consists  of  numerically 
integrating  the  equations  of  motion  to  solve  for  the  trajectory  at 
discrete  time  points.  The  classical  statistical  mechanics  relates 
thermodynamic  properties  to  the  average  value  of  functions  of  the 
phase  coordinates  of  the  atoms.  The  average  is  usually  taken  over 
all  elements  of  the  ensemble  which  corresponds  to  a  given  macrostate. 

The  present  work  is  a  study  of  the  static  properties  of  fracture 
in  crystalline  solids.  Using  CMD  simulation  results,  a  number  of 
continuum  theory  predictions  are  tested  at  an  atomic  level,  and  the 
effects  of  non-linear  relaxation  observed. 

It  is  well  known  that  to  find  a  general  solution  to  the  static 
configuration  around  defects  of  irregular  geometries  can  be  a  dif¬ 
ficult  task  even  in  the  approximation  of  a  linear  theory.  Moreover, 
there  are  cases  where  the  system  behavior  is  determined  by  a  region 
that  is  only  a  few  interatomic  spacings  in  extent;  then  non-linear 
atomistic  processes  can  not  be  ignored.  In  such  cases  CMD  can  be 
an  extraordinarily  useful  technique  for  studying  the  static  configura¬ 
tion  of  the  defect  as  well  as  its  dynamic  properties. 
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The  first  studies  of  fracture  in  a  discrete  lattice  were  car¬ 
ried  out  with  very  simple  models  in  one  or  two  dimensions  and  highly 
idealized  interatomic  force  laws.^  These  simplified  models,  how¬ 
ever,  contributed  greatly  to  the  understanding  of  the  macroscopic 
properties  of  fracture,  such  as  brittle  versus  ductile  behavior,  and 
determination  of  elastic  constants,  ideal  maximum  cohesive  and  shear 
stress,  etc.,  in  terms  of  microscopic  physical  variables,  such  as 
lattice  structure,  interatomic  force  law,  etc. 

Subsequent  studies  using  computer  simulation  were  concerned 
with  systems  of  larger  number  of  particles  and  more  realistic  inter¬ 
atomic  force  laws,  such  as  the  calculations  for  iron  by  Kanninen  and 
Gehlen^  and  for  silicon  by  Sinclair  and  Lawn.  ^ 

The  main  difficulty  with  these  computer  simulations  has  been 
the  specifications  of  the  external  boundary  conditions.  Periodic 
boundary  conditions  is  a  simple  and  accurate  way  of  simulating  the 
properties  of  an  infinite  homogeneous  system.  However,  when  the 
system  is  acting  under  the  effect  of  an  external  force,  a  different 
procedure  must  be  used.  One  alternative  is  to  apply  the  external 
forces  directly  over  the  free  surfaces  of  the  system.  The  other 
consists  of  matching  the  solution  at  the  boundaries  to  a  continuum 
elasticity  solution.  Of  the  two  alternatives,  only  the  second  has  a 
practical  application  in  the  study  of  a  macroscopic  crack  under  the 
action  of  an  exterior  stress  applied  far  from  the  crack  and  it  alone 
will  be  used  in  our  simulations. 

Dynamics  of  crack  propagation  has  been  studied  by  computer 
molecular  dynamics  by  Sander, ^  by  Weiner  and  Pear ^  and  by  Ashurst 
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and  Hoover.  All  these  studies  show  the  influence  of  atomic  struc¬ 
ture  and  interatomic  force  laws  on  the  dynamics  of  crack  propagation. 
These  simulations  used  free  boundaries  and  the  stresses  were  applied 
directly  over  the  exterior  atoms,  close  to  the  crack  tip.  It  was 
found  that  the  results  differ  substantially  from  predictions  given  by 
macroscopic  theories,  and  it  seems  appropriate  to  consider  more 
realistic  border  conditions  if  one  wants  to  simulate  the  behavior  of 
cracks  in  real  materials  or  establish  a  direct  comparison  between 
continuum  and  atomistic  results. 

Our  study  of  static  properties  of  fracture  by  CM)  consists  of 
simulating  the  crack  tip  behavior  in  a  two-dimensional  atomic  system, 
introduced  in  an  infinite  continuum  medium  under  the  action  of  an 
external  stress  (mode- I  fracture) . 

The  static  equilibrium  configuration  of  the  system  is  obtained 
by  applying  a  frictional  damping  to  the  equations  of  motion.  The 
general  techniques  used  in  this  work  jointly  with  the  standard  methods 
of  CMD  are  summarized  in  Chapter  Two,  where  a  description  of  the 
problems  associated  with  the  different  types  of  boundary  conditions 
is  also  included. 

Chapter  Three  summarizes  several  important  concepts  in  con¬ 
tinuum  theory  of  fracture,  namely,  the  Griffith's  theory  of  brittle 
materials,  and  the  criteria  which  have  been  recently  developed  to 
predict  brittle  or  ductile  behavior  in  crystalline  materials.  The 
last  section  of  this  chapter  describes  previous  research  done  in 


this  area  by  CMD. 
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The  description  of  the  technique  used  in  this  work,  flexible 
boundary,  is  done  in  Chapter  Four,  jointly  with  the  method  to  de¬ 
termine  the  elastic  constants,  ideal  stresses  and  surface  energy  of 
a  material  whose  atoms  interact  through  a  known  pair  potential 
function.  These  physical  constants  are  important  not  only  for  cal¬ 
culating  the  displacement  field  according  to  linear  continuum  solu¬ 
tions,  but  also  for' predicting  brittle  or  ductile  behavior  of  the 
system.  In  the  last  part  of  the  chapter  the  critical  stress  to 
propagate  a  crack  in  a  given  material  is  determined  by  direct  appli¬ 
cation  of  the  flexible  boundary  technique  and  these  results  are 
compared  with  continuum  linear  elasticity  theory. 

Plasticity  around  the  crack  tip  has  been  studied  in  Chapter 
Five.  Two  criteria  to  determine  brittle  or  ductile  behavior  of  a 
given  material,  one  due  to  Kelly,  Tyson  and  Cottrell  (KTC)  and 
another  due  to  Rice  and  Thomson  (RT) ,  are  studied  and  their  predic¬ 
tions  compared  to  CMD  simulations. 

Emission  of  dislocations  and  several  microscopic  processes  of 
plastic  deformation  have  been  observed  in  the  simulation  results. 

The  strain  and  rotation  fields  around  the  crack  tip  for  several  poten¬ 
tial  functions  have  been  determined.  Finally,  a  summary  with  conclu¬ 
sions,  discussion  of  results,  and  suggestions  of  further  work  are 


given  in  Chapter  Six. 
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2.1  Equations  of  Motion 

In  CMD  calculations,  a  three-dimensional  system  of  N  particles 
is  treated  by  setting  up  3  N  classical  equations  of  motion  which  are 
coupled  through  an  assumed  two-body  interaction  potential.  This 
set  of  3  N  simultaneous  differential  equations  is  then  integrated 
numerically  on  a  computer  to  give  the  spatial  trajectories  and 
velocities  of  all  the  particles  as  a  function  of  time. 

In  the  present  simulation,  a  classical  system  of  N  point  par¬ 
ticles  of  mass  m  is  assumed  to  obey  Newtonian  mechanics,  in  which 
case  the  equations  of  motion  are  given  by 

-  J-  h  (  K  ^  r»  ]  +■  j  N  (*‘0 

cU1  '  ™  ^ 


where  is  the  portion  coordinate  of  the  ith  particle,  and  F^ 
are  the  vector  forces  on  the  ith  particle  due  to  other  particles  and 
exterior  force  respectively.  For  central,  conservative,  pairwise 
additive  potentials,  the  pair  force  F_  between  two  particles  i  and 
j  of  the  system  is  the  gradient  of  j>  (  'Qj)  with  respect  to  ,  i.e.. 


(2.-2.) 


where  ^  (  'Qj  )  is  the  pairwise  potential  between  the  particles  i  and  j. 
It  is  a  function  of  the  pair  separation  vector  between  these  two 


particles. 
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2.2  Potential  Function 

Empirical  and  semi-empirical  potential  functions  can  be  con¬ 
structed  in  which  adjustable  parameters  are  determined  by  matching 
certain  calculated  properties  of  the  system  to  experimental  results. 
Two  fundamental  assumptions  are  almost  universally  employed:  central 
forces  and  pairwise  additivity.  Pairwise  additivity  means  that  the 
potential  energy  between  two  atoms  is  unaffected  by  the  presence  of 
other  nearby  atoms.  The  central  force  assumption  means  that  the 
force  between  atoms  is  directed  along  the  line  joining  the  centers 
of  mass  and  so  is  a  function  only  of  the  atomic  separation.  The 
next  discussion  goes  into  commonly  known  empirical  and  theoretical 
potential  functions  which  have  been  used  in  our  computer  simulations. 

Some  interatomic  potentials  are  characterized  by  a  single 
functional  form  which  varies  continuously  over  the  entire  range  of 
interatomic  spacing  and  has  short  range  repulsive  and  long  range  at¬ 
tractive  forces.  A  potential  of  this  kind  was  employed  by  Lennard- 
Jones^  in  which  the  attractive  and  repulsive  parts  vary  as  the 
inverse  of  the  12  and  6  power,  respectively.  This  potential  has 
often  been  used  to  represent  the  effective  pair  interactions  in  rare 
gases  and  in  simple  non-polar  molecules  as  and  CH^: 

This  potential  gives  qualitative  agreement  with  rare  gas  properties 
in  condensed  phases  and  is  computationally  convenient  because  the 
entire  potential  is  expressed  in  pairwise  additive  two-body  terms. 
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(2) 

A  form  extensively  used  for  metals  is  the  Morse  potential 


r  -Lot 
<p(r)z  Dj  H 


Cv--ir„)  .«(*-.  r.) 

-10, 


where  «  and  D  are  constants  with  dimensions  of  reciprocal  distance 
and  energy,  respectively,  and  ^  is  the  equilibrium  distance  cor¬ 
responding  to  the  minimum  of  the  potential.  The  Morse  potential 
approaches  zero  exponentially  at  large  distance,  in  accord  with  the 
Thomas-Fermi  model  of  electronic  screening.  It  also  yields  an  ex¬ 
actly  soluble  Shrodinger  equation  for  the  vibrational  eigenvalue 

problem  of  diatomic  molecules  and  the  eigenvalue  distribution  is  in 

(2) 

good  agreement  with  experimental  spectroscopic  measurement. 

The  most  important  characteristics  of  this  potential  pertaining 

to  our  CMD  simulation  of  fracture  is  that  the  three  parameters  D, 

d  and  rD  of  the  potential  for  several  f.c.c.  and  b.c.c.  metals  have 

been  determined  using  experimental  values  for  the  energy  of  vapori- 

(3) 

zation,  the  lattice  constant,  and  the  compressibility. 

The  estimates  of  the  validity  of  the  Morse  potential  in  metals 
given  in  the  preceding  paragraphs  apply  to  a  perfect  crystal.  An 
additional  complication  is  introduced  if  defects  are  present.  The 
Morse  potential  constants  computed  from  the  energy  of  vaporization, 
the  lattice  constant,  and  the  compressibility  refer  to  a  perfect 
lattice  and  reflect  its  electronic  distribution.  The  presence  of  a 
defect,  however,  alters  the  electron  distribution,  and  it  is  diffi¬ 
cult  to  estimate  how  this  altered  distribution  would  affect  the  atoms 

(3) 

in  the  vicinity  of  the  defect. 
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2.3  Numerical  Integration  of  Equation  of  Motion 

The  algorithm  used  to  move  the  particles  forward  in  time  by  an 
amount  At  is  given  by 


rC  t  +  ab)  5  rit)  +  vCtJat  +  JL  [  4  a.lfc)  -a.((r-at)3 

Cl-5) 


\z(t  +  afc)r  v (t)  -t-  +  +  ~  0>(t~Ak)J 

6 


where  r(t),  v(t)  and  a(t)  are  the  position,  velocity  and  acceleration 
of  a  particle  at  time  t.  This  algorithm,  used  by  P.  Shofield  in  CMD 

(4) 

simulation  of  liquids  is  very  strongly  conserving  and  allows  a 
relatively  large  time  step  to  be  used. 

We  have  checked  the  accuracy  of  this  algorithm  by  comparing  it 
with  the  exact  solution  for  the  harmonic  oscillator.  Both  solutions, 
the  exact  and  the  approximation,  have  been  expanded  in  powers  of  t. 

(a)  Exact  solution  corresponding  to  an  harmonic  oscillator 
XCb  -  xC o)  v  Y±£J  -  V(0)t  i  Cwt) 

VAl 

Substituting  t  by  t  +  at  in  the  above  equations  and  expanding 
Si.n  0*4  + at))  OLoiol  Ct+afcjJ 

x(fc  +  ah  -  xfbWvffcl&b  f  afr  Q-Ct)  -  &b  +  V,  jiVmjt 

2-  3  I  &  ^  i  (Z.-6) 

v(.t  +&t )  vtb)  +  cvU)at  +  v0c*ioofc  y,  va 

l'  3!  ~T~I 
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(b)  Algorithm  (2-5) ,  expanded  in  powers  of  a  t 


i  ! 


)  2 


v(fc  +aI)-.  y(t)  4-  a.[t)*k  +  v0c*ut('^_£jL)  +  v«  ^ 

2-  '■  '  II 


(2-7J 

i 

i... 


Equation  (2-7)  shows  that  the  algorithm  used  by  Schoefield  is  cor- 

3  2 

rect  to  order  (At)  and  (At)  for  the  positions  and  velocities,  re¬ 
spectively,  when  it  is  tested  against  the  harmonic  oscillator 
solution. 

This  algorithm  has  the  advantage  that  it  needs  less  terms  to 
calculate  r(t  +  a  t)  and  v  (t  +it)  in  expression  (2-5)  than  the  ordinary 
central  difference  method  for  the  same  order  of  accuracy  and  the  dis¬ 
advantage  that  the  evaluation  at  time  t  +4  t  is  carried  out  using  the 
values  at  t  and  t  -  At,  therefore  a  different  procedure  is  needed  to 
initialize  the  numerical  integration  at  time  t  =  0. 

2.4  Accuracy  of  Numerical  Integration;  Conservation  of  Energy  and 

Reversibility  of  Trajectories 

The  time  step  size  has  to  be  chosen  small  enough  that  the  inte¬ 
gration  procedure  generates  a  stable  solution  to  the  equations  of 
motion.  If  the  time  step  size  is  too  large,  then  the  equations  of 
motion  will  fail  to  conserve  total  energy  and  will  result  in  an  un¬ 
stable  phase  space  trajectory.  On  the  other  hand,  it  also  has  to  be 
chosen  as  large  as  possible  because  the  computational  time  per  time 
step  does  not  depend  on  the  time  step  size,  and  the  net  cost  of  a 
simulation  over  a  fixed  duration  will  vary  inversely  as  the  time  step 
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size.  The  numerical  integration  should  also  span  an  oscillation 
period  in  a  number  of  time  steps,  i.e.,-^10,  such  that  the  net  change 
during  any  single  time  step  is  not  too  large.  The  accuracy  of  the 
numerical  integration  can  be  measured  by  the  precision  to  which  the 
total  energy  is  conserved.  A  series  of  simulations  with  different 
time  step  sizes  at  different  temperatures  was  carried  out  on  a  per¬ 
fect  lattice  composed  of  64  particles  to  determine  the  maximum  time 
step  size  that  gave  reasonable  conservation  in  the  total  energy.  Time 
step  sizes  as  large  as  1/10  of  an  oscillation  period  gave  absolute 
deviations  of  about  0.05%  per  time  step,  and  were  considered  accept¬ 
able.  However,  more  conservative  time  step  sizes,  about  1/50  of  one 
oscillation  period,  were  used. 

The  accuracy  of  the  numerical  integration  can  also  be  measured 
by  the  accuracy  of  the  reversibility  of  trajectory.  The  trajectory 
can  be  retraced  simply  by  reversing  velocities  and  other  odd  time 
derivations  at  the  time  step  of  retracing  or  well  by  using  a  negative 
value  f or  a t  at  every  time  step.  The  accumulated  error  in  the  total 
energy  during  a  computation  period  of  one  thousand  time  steps  (500 
time  steps  before  reversing  trajectories),  at  a  temperature  of  60°K, 
was  less  than  2%. 

2.5  Potential  Range  Cut-off 

The  number  of  neighbors  that  are  included  in  a  CMD  simulation 
is  an  important  parameter  because  the  calculation  time  scales  linearly 
with  the  total  number  of  interacting  neighbors.  The  total  number  de¬ 
pends  on  the  potential  cut-off  or  the  distance  at  which  the  force  be¬ 
comes  negligible  by  comparison  with  the  force  between  nearest  neighbors. 
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However,  the  elimination  of  direct  interaction  between  a  pair  of 
particles  separated  by  other  interacting  particles  does  not  preclude 
indirect  communication  between  those  particles.  Each  particle  in 
the  system  is  dynamically  coupled  to  every  other  particle  in  the 
system  either  directly  or  indirectly. 

In  the  simulation  of  fracture  the  cut-off  range  plays  an  imr 
portant  role  because  the  propagation  of  the  crack  is  produced  by 
bond  rupture  in  tension  when  they  reach  the  cut-off  range,  or  by 
bond  rupture  in  shear  in  which  case  the  crack  is  blunted.  The  surface 
energy  is  another  parameter  that  can  change  significantly  if  an  ap¬ 
propriate  cut-off  range  is  not  chosen.  We  have  chosen  in  our  simu¬ 
lations  only  first  nearest  neighbor  interaction  with  the  potential 
function  truncated  at  a  distance  where  the  force  can  be  neglected. 
Besides,  several  tests  were  done  using  second  nearest  neighbors  and 
the  results  did  not  change  significantly. 

2.6  Boundary  Conditions 

The  behavior  of  an  infinite  system  is  approximated  by  describ¬ 
ing  the  motion  of  N  particles  in  a  finite  computational  cell  of 
volume  V  with  periodic  boundary  conditions. 

If  one  wants  to  use  this  procedure  to  study  the  behavior  of 
defects,  the  system  must  be  large  enough  such  that  the  images  do  not 
interact  with  the  original  defect.  An  additional  difficulty  arises 
when  the  defect  is  subject  to  a  certain  state  of  stress.  Then,  the 
position  of  the  particles  must  be  generated  according  to  some  known 
theory  specifying  that  state  before  periodic  boundary  conditions  can 


be  applied. 
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To  overcome  these  difficulties  in  the  simulation  of  fracture, 
two  general  procedures  have  been  implemented  during  the  last  few  years. 
One  consists  of  generating  the  state  of  stress  by  applying  forces 
directly  over  the  free  boundaries  of  the  system  with  the  defect.  The 
disadvantage  of  using  this  procedure  is  that  it  is  very  expensive  to 
work  with  a  system  large  enough  to  avoid  interactions  between  bound¬ 
aries  and  the  defect.  The  other  method  arranges  the  atoms  at  the 
boundaries  according  to  the  solution  given  by  the  linear  elasticity 
theory  and  studies  the  behavior  within  the  system,  where  linear  elas¬ 
ticity  theory  is  in  general  not  valid,  using  computer  molecular 
dynamics . 

The  last  procedure  which  is  the  one  used  in  this  work  has  the 
great  advantage  of  reproducing  exactly  the  macroscopic  conditions 
with  a  relatively  small  number  of  particles  letting  exact  calculations 
of  the  defect  structure.  This  procedure  is  called  fixed  boundary 
condition. 

A  further  improvement  appears  when  it  is  found  that  after  the 
relaxation  of  the  atomic  non-linear  region  using  CMD,  small  residual 
forces  arise  at  the  boundaries  due  to  the  difficulties  of  the  con¬ 
tinuum  linear  solution  applied  at  the  boundaries  to  accomodate  the 
atomic  region  after  the  relaxation.  In  this  case,  what  has  been 
recently  applied  is  a  continuum  linear  Green's  function,  which  gives 
the  displacement  field  due  to  a  unit  force  at  some  distance  from  the 
crack  necessary  to  cancel  the  small  remaining  forces. 

This  method  called  flexible  boundary  has  been  applied  in  this 
work  to  determine  crack  tip  configuration  and  the  critical  stress 
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intensity  factor  of  brittle  materials.  In  a  few  cases  where  it  was 
found  a  large  plastic  region,  formed  by  emission  of  dislocations  from 
the  crack  tip,  the  first  steps  of  the  technique,  which  constitutes 
the  fixed  boundary  method,  was  applied  to  a  large  system  to  avoid  the 
interaction  of  the  extending  plastic  region  with  the  boundaries. 

2.7  Determination  of  Equilibrium  Configuration 

Several  different  kinds  of  computer  experiment  methods  can  be 
used  to  obtain  the  defect  configuration.  They  are  classified,  ac¬ 
cording  to  Beeler, ^  in  five  types:  dynamical,  frictional  damping, 
variational,  Monte-Carlo  and  lattice-static  methods.  A  pairwise 
atomic  interaction  model  is  assumed  for  computing  structure- dependent 
forces  and  energy  contributions.  Structure  independent  forces  and 
their  associated  energy  contribution  usually  are  accounted  for  by 
specification  of  appropriate  boundary  forces  and  the  atomic  volume. 
Detailed  descriptions  of  the  first  four  methods  and  their  applications 
are  given  in  Ref.  6.  The  lattice-static  method  is  explained  in 
Refs.  7-9. 

In  the  dynamical  simulation  of  a  crystal,  each  atom  vibrates 
continuously  about  its  static  equilibrium  position.  The  array  of 
atomic  positions  can  be  obtained  by  computing  the  time-average  portion 
of  each  atom.  In  instances  where  only  the  time-average  portions  are 
desired,  the  dynamical  method  can  be  amended  so  that  the  atom  motion 
is  progressively  damped  as  the  atoms  approach  their  static  equilibrium 
portions . 

There  are  several  procedures  for  reaching  the  static  equilibrium 
configuration.  Gibson  et  al.  introduced  a  frictional  method  for 
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static  equilibrium  configuration  calculations  wherein  the  velocity 
of  each  atom  is  set  to  zero  each  time  the  kinetic  energy  of  the  en¬ 
tire  system  reaches  a  maximum.  This  procedure  is  suggested  by  the 
circumstance  that  the  velocity  magnitude  of  a  mass  point  in  simple 
harmonic  motion  is  maximum  at  the  time  it  is  passing  through  the 
static  equilibrium  position  (zero  net  force  position).  Several  vari¬ 
ations  of  this  basic  idea  have  evolved.  As  a  mass  point,  in  simple 
harmonic  motion,  approaches  the  static  equilibrium  position,  the  dot 
product  of  its  velocity  and  the  restoring  force  in  a  monotomically 
decreasing  positive  quantity  which  vanishes  at  the  equilibrium  posi¬ 
tion  and  then  becomes  negative.  On  this  basis,  Evans  and  Beeler 
introduced  the  idea  of  setting  the  velocity  of  any  atom  to  zero 
whenever  the  velocity-force  dot  product  becomes  negative,  in  a  static 
equilibrium  calculation.  This  individual  atom  criterion  for  damping 
leads  to  a  significantly  faster  convergence  than  the  maximum  total 
kinetic  energy  criterion. 

The  procedure  used  in  this  work  consists  of  introducing  a 
frictional  force  in  the  equations  of  motion  of  each  particle  after 
the  system  has  released  freely  for  several  oscillations  of  the  total 
kinetic  energy: 

Fc  (  r*  ,  , . . ..  tr„)  -  *  <LH  (  l- 

d  «*■  fc 

The  damping  constant  has  to  be  chosen  small  enough  that  the  system 
is  not  so  overdamped  that  it  is  quenched,  and  as  large  as  possible 
so  the  system  will  reach  the  minimum  internal  energy  configuration 
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in  a  reasonable  period  of  time.  Usually  the  damping  constant  is 
chosen  to  be  around  half  the  time  average  value  of  the  critical  damp¬ 
ing  constant  estimated  by  Hooke's  Law.  We  did  several  tests  with 
different  damping  constants  and  found  that  for  values  close  to  the 
critical  damping,  the  average  time  for  the  system  to  reach  the  minimum 
internal  energy  configuration  was  about  300  time  steps.  The  total 
procedure,  therefore,  consists  of  several  periods  of  free  relaxation 
alternating  with  periods  with  frictional  damping  until  the  system 
temperature  has  decreased  to  a  few  degrees  Kelvin.  The  final  con¬ 
figuration  corresponds  to  a  stable  state  of  minimum  potential  energy 
and  further  damping  to  still  lower  temperatures  does  not  lead  to  any 
significant  different  configuration  of  the  relaxed  atomic  positions. 


Chapter  Three 


Computer  Molecular  Dynamic  Simulation  of  Fracture 


3.1  Griffith  Theory  of  Brittle  Materials 

3.2  Brittleness  and  Ductility  of  Crystalline 

Materials 

3.3  Computer  Molecular  Dynamics  of  Fracture  — 

Previous  Works 
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The  following  sections  are  devoted  to  summarize  some  of  the 
general  concepts  and  theories  of  macroscopic  fracture  mechanics  which 
will  be  examined,  from  an  atomistic  point  of  view,  in  the  present  work. 
A  description  of  the  different  models  and  recent  results  in  atomic 
simulation  of  fracture  is  given  in  Section  3.3  of  this  chapter. 

3.1  Griffith  Theory  of  Brittle  Materials 

Griffith's  idea  was  to  set  up  a  model  for  a  crack  system  in 

(1  2) 

terms  of  a  reversible  thermodynamical  process.  5  The  important 
elements  of  the  system  are  defined  in  Fig.  3.1:  an  elastic  body  B 
containing  an  internal  crack  S  of  length  2C  is  subjected  to  loads 
applied  at  the  outer  boundary  L. 

The  first  step  in  the  treatment  was  to  find  an  expression  for 
the  total  energy  of  the  system.  To  do  this  the  individual  energy 
terms  which  change  as  a  result  of  crack  formation  are  considered. 
First,  it  can  be  expected  in  general  that  the  outer  boundary  of  the 
cracked  body  will  undergo  some  displacement,  such  that  the  applied 
load  does  an  amount  of  work  W  .  (For  a  truly  reversible  system  an 
increase  in  this  work  can  be  identified  with  a  corresponding  decrease 
in  the  potential  energy  of  the  loading  system.)  Second,  the  strain 
energy  V,,  stored  in  the  elastic  medium  must  be  sensitive  to  any  vari- 
ations  in  the  system  geometry.  Third,  the  mere  act  of  creating  new 
crack  surfaces  requires  the  expenditure  of  free  surface  energy  Vg. 

For  a  static  crack  system  the  total  energy  is  the  sum  of  threse 
three  terms: 

(3-/) 


V=  (-K  -Vf)  +vs 


Figure  3.1 


Static  plane  crack  system.  L  applied 
loading,  S  crack  surface. 


C  £<*<a.ck  Jle/v^tk) 


Figure  3.2 


Energetics  of  Griffith  crack  in 
uniform  tension. 
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Since  forces  transmitted  to  the  crack  region  are  determined  by  the 
loading  system  and  elastic  medium,  it  is  convenient  to  refer  to  the 
bracket  term  in  Eq.  (3-1)  as  the  mechanical  energy  of  the  system. 

Of  course,  if  we  were  to  be  concerned  with  a  dynamic  crack  system, 
we  would  have  to  add  a  kinetic  energy  term  to  Eq.  (3-1). 

Thermodynamic  equilibrium  is  then  attained  by  balancing  the 
mechanical  and  surface  energy  terms  over  a  virtual  crack  extension 
.  The  mechanical  energy  must  decrease  as  the  crack  extends.  On 
the  other  hand,  the  surface  energy  term  must  increase  with  crack 
extension,  since  the  cohesive  forces  of  molecular  attraction  across 
Sc  must  be  overcome  during  the  creation  of  the  new  fracture  surfaces. 
Thus,  the  bracket  term  in  Eq.  (3-1)  favors  crack  extension,  while  the 
second  opposes  it.  This  is  the  Griffith  energy  balance  concept,  a 
formal  statement  of  which  is  given  by  the  extended  equilibrium 
requirement 


ciV  . 

<j-  C 


(3-Zj 


Griffith  attempted  to  confirm  his  theory  by  applying  it  to  a 

real  situation.  He  needed  a  model  for  a  crack  in  order  to  calculate 

the  energy  terms  in  Eq.  (3-1).  For  this  he  took  advantage  of  the 
(3) 

Inglis  analysis,  considering  the  case  of  a  narrow  elliptical  crack 
in  a  remote,  uniform  tensile  stress  field.  Then,  for  experimental 
verification  he  had  to  find  a  well  behaved,  model  material,  isotropic 
and  closely  obeying  Hooke's  law  at  all  stress  levels  prior  to  fracture. 
Glass  was  selected  as  the  most  easily  accessible  material  satisfying 


these  demands. 
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In  evaluating  the  mechanical  energy  of  his  model,  Griffith  used 
a  result  from  linear  elasticity  theory,  namely  that  for  any  body  under 
constant  applied  stress  during  crack  formation, 

\fl/  -  V£  (  ceY^-to^t  J  (  3  •  3 J 

From  the  Inglis  solution  of  the  stress  and  strain  fields,  the 

strain  energy  density  is  readily  computed  for  each  volume  element 

about  the  crack.  Integrating  over  dimensions  large  compared  with  the 

length  of  the  crack  then  gives,  for  unit  width  of  the  crack  (measured 

(3  4) 

along  the  crack  front  ’ 

V  -  7T  C  j  E  piasrJ-  aLmlu  J  (3 -4  a.} 

for  a  thin  plate,  or 

V  z  7 T  (i  -  Cl  $~Ll  /  £  ( ft'iO**)  (i  -  4  b) 

for  a  thick  plate.  Here  XTU  is  the  applied  tension  normal  to  the 
crack  plane,  E  is  Young's  modulus,  V  is  Poisson's  ratio,  and  2C  is 
the  length  of  the  crack.  The  application  of  additional  loading 
parallel  to  the  crack  plane  has  negligible  effect  on  the  strain  energy 
term  in  Eq.  (3-4).  For  the  surface  energy  of  the  crack  system, 
Griffith  wrote  for  unit  width  of  crack 


v}  --  ^  <  r 


C  i- 5) 
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with  ^  the  free  surface  energy  per  unit  area.  The  total  system  energy 
Eq.  (3-1)  thus  becomes  for  the  case  of  plane  stress 

V:  -TT  CLTlL  j  E  f  ^  C  f 

The  Griffith  equilibrium  condition,  Eq.  (3-2),  may  now  be  applied  to 
Eq.  (3-6);  this  gives  as  a  critical  condition  for  fracture 

(T  r  E  ^  /  Tr  d  )  (3~7) 

for  constant  stress,  plane  stress  conditions  (Fig.  3.2). 

In  the  above  case,  the  analysis  has  been  applied  to  a  crack  in 
a  brittle  elastic  solid  where  its  movement  does  not  involve  plastic 
deformation.  It  is  significant  that  the  same  approach  can  also  be 
applied  to  problems  where  the  movement  of  the  crack  involves  plastic 
deformation  with  the  assumption  that  the  plastic  region  around  the 
crack  tip  is  negligibly  small  in  comparison  with  the  outer  zone 
(Irving  approach) .  In  this  case  we  may  substitute  by  ,  where 
is  the  plastic  work  done  in  forming  a  square  centimeter  of  surface. 
Under  these  circumstances  the  ^  surface,  defined  in  terms  of  the 
energy  to  break  atomic  bonds  must  include  the  energy  due  to  all  the 
mechanism  of  plastic  relaxation  produced  during  the  crack  extension. 
3.2  Brittleness  and  Ductility  of  Crystalline  Materials 

There  is  extensive  literature  covering  all  the  continuum  and 
atomistic  aspects  of  crack  tip  plasticity.  In  this  section  a  model 
of  dislocation  nucleation  at  the  tips  of  cracks  and  the  influence  of 
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such  defects  to  predict  brittle  or  ductile  behavior  as  described  by 
Kelly,  Tyson  and  Cottrell  (KTC) ^  and  by  Rice  and  Thomson  (RT) ^ 
is  discussed. 

A  criterion  for  brittle  fracture  in  crystals  can  be  established 

in  terms  of  the  spontaneous  emission  of  dislocations  from  a  sharp 

crack.  The  stress  field  near  the  crack  tip  in  a  linearly  elastic 

-’4  , 

medium  is  of  the  form  j  K  (iff'C)  l  © )  where  K  is  the  stress 
intensity  factor  and  (r,8)  are  polar  coordinates  with  the  origin  at 
the  crack  tip.  At  the  crack  tip  itself  where  Y~-»  0  a  non-linear 
treatment  is  required  as  the  interatomic  bonds  are  stretched  beyond 
the  region  of  harmonic  behavior.  As  the  stress  intensity  increases, 
lattice  failure  will  occur  at  the  crack  tip  in  the  non-linear  region. 
Failure  can  occur  by  bond  rupture  in  either  tension  or  shear,  which 
will  determine  whether  the  behavior  is  inherently  brittle  or  ductile. 

Two  treatments  have  been  presented  in  the  literature  to  predict 
the  failure  mode  from  known  material  properties.  Kelly,  Tyson  and 
Cottrell  proposed  the  criterion  that  an  ideal  (defect-free)  solid 
may  sustain  a  fully  brittle  crack  only  if  the  theoretical  strength 
in  tension  is  exceeded  before  the  theoretical  strength  in  shear  within 
the  local  field  of  the  tip.  The  criterion  makes  no  comment  on  the 
nature  of  any  shear  deformation  that  might  happen.  One  only  needs  to 
consider  the  theoretical  strength  calculations  in  relation  to  the 
stress  distribution  at  the  crack  tip,  coming  back  to  the  equations 
of  linear  elasticity  for  an  evaluation  of  this  field.  The  final  ex¬ 
pression  proposed  by  KTC  was  that  brittle  fracture  would  be  observed 


where  the  subscripts  "ideal"  and  "max" 
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refer  to  the  ideal  properties  of  a  perfect  lattice  and  to  the  maximum 
values  attained  at  a  crack  tip  respectively. 

Rice  and  Thomson on  the  other  hand,  have  argued  that  a  neces¬ 
sary  criteron  for  brittle  fracture  is  stability  against  the  emission 
of  dislocations  from  the  crack  tip.  They  have  treated  this  nucleation 
process  within  the  approximation  of  linear  elasticity  and  the 
model  of  a  dislocation  core,  the  specific  question  addressed  by  these 
workers  is  whether  the  shear  deformation  is  sufficient,  not  merely  to 
nucleate  a  dislocation  at  the  crack  tip,  but  also  to  propagate  it 
within  the  stress  field  into  the  surrounding  crystal.  They  considered 
first  the  balance  between  the  following  crack-dislocation  interaction 
forces:  (a)  force  on  dislocation  arising  from  stress  field  of  crack, 
(b)  surface  tension  force  caused  by  creating  surface  ledges  at 
(blunted)  crack,  (c)  image  force  of  dislocation  in  the  free  surface 
of  crack.  The  first  term  repels  the  dislocation  from  the  crack  tip, 
while  the  remaining  two  attract  it.  One  then  compares  the  cores  of 
the  dislocation  versus  the  distance  in  which  the  dislocation  starts 
to  be  repelled  from  the  crack  tip. 

Rice  and  Thomson  thus  conclude  that  covalent  and  ionic  solids, 
along  with  h.c.p.  metals,  are  stable  against  dislocation  emission, 
while  f.c.c.  metals  are  unstable;  b.c.c  metals  comprise  an  inter¬ 
mediate  case.  The  KTC  criterion  probably  leads  to  an  underestimate 
of  the  brittle  tendencies  of  the  solids  because  of  not  considering 
the  subsequent  propagation  of  newly  generated  shear  deformation 
(dislocation) . 
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Both  treatments  require  approximations  in  their  quantitative 

development.  The  stress  analysis  of  linear  elasticity  is  used  to 

evaluate(GJht/i/\ /57ch£-14Vg)  around  the  crack  tip  in  the  KTC  criterion  and 

the  forces  on  dislocations  near  the  crack  tip  in  the  RT  model. 

Better  models  of  the  non-linear  and  atomistic  behavior  of  materials 

are  also  required  to  evaluate  )  and  dislocation  core 

idlfll 

structure  with  greater  accuracy. 

3.3  Computer  Molecular  Dynamics  of  Fracture  -  Previous  Works 

During  the  last  years  a  number  of  studies  of  fracture  in  crystal¬ 
line  solids  using  atomistic  models  have  appeared.  One  can  classify 
them  in  two  categories.  The  first  type  is  concerned  with  a  specific 
material  or  group  of  materials,  and  its  purpose  is  to  be  as  realistic 
as  possible  with  respect  to  both  crystal  geometry  (dimensionality  of 
the  system,  boundary  conditions,  etc.)  and  interatomic  force  laws. 

These  models  must  yield  information  regarding  atomic  processes  in  the 
given  class  of  materials.  Examples  are  provided  by  the  work  of 
Kanninen  and  Gehlen^’^’^  and  Sinclair. 

Gehlen  and  Kanninen  worked  on  oc-iron  in  a  three-dimensional 
lattice  on  which  the  atoms  interact  through  a  central  pair  potential 
constructed  by  Johnson.  This  potential  is  based  upon  two-body  inter¬ 
actions  which  extend  out  to  first  and  second  nearest  neighbors.  An 
empirical  form  is  used  with  the  constants  selected  to  match  the  experi¬ 
mentally  determined  elastic  constants  of  the  material.  The  atoms  are 
initially  put  into  a  configuration  approximating  the  defect  being 
simulated.  The  boundary  atoms  are  held  fixed  in  the  positions  given 
by  the  linear  elastic  continuum  solution  (fixed  boundary  condition) . 
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The  final  position  of  atomic  equilibrium  is  evaluated  by  letting  the 
atoms  inside  the  system  move  freely  according  to  the  interatomic 
forces. 

Based  on  this  model  Gehlen  and  Kanninen  tried  to  check  the 
validity  of  Griffith's  criterion.  All  the  macroscopic  variables, 
elastic  constants  and  surface  energy  were  specified  by  the  potential 
function  and  the  critical  stress  was  determined  by  calculating  the 
condition  of  maximum  elongation  for  the  crack  tip  bond.  The  estimated 
values  for  the  critical  stress  intensity  factor  were  found  to  be 
greater  than  those  necessary  to  produce  an  exact  correspondence  with 
Griffith's  equation.  A  possible  reason  for  this  discrepancy,  as 
noted  by  the  authors,  is  that  the  boundaries  were  fixed.  This  inflexi¬ 
bility  could  produce  high  residual  forces  between  the  atomic  and  con¬ 
tinuum  regions,  especially  if  the  atomic  boundaries  are  close  to  the 
crack  tip.  To  avoid  these  difficulties  the  authors  used  a  new  pro¬ 
cedure  called  flexible  boundary  to  attenuate  the  effects  of  the  lack 
of  adjustment  between  the  two  regions.  In  the  present  work  an  im¬ 
provement  of  this  method  is  applied  to  crack  tips  and  it  will  be  dis¬ 
cussed  in  Chapter  Four. 

The  second  type  of  atomistic  models  is  highly  idealized  from 
the  viewpoint  of  crystal  geometry  (two  dimensions,  free  boundaries, 
etc.)  and  interatomic  force  laws.  They  cannot,  therefore,  represent 
directly  any  particular  real  material.  Nevertheless,  they  do  serve 
to  provide  insight  into  general  characteristics  of  static  and  dynamic 
propagation  of  crack  with  a  minimum  of  computational  difficulties. 
Examples  of  idealized  atomistic  models  include  the  works  of  Thomson, 
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Hsieh  and  Rana^11,12,13'*  in  which  a  new  phenomenon  called  "lattice 
trapping",  the  crack  analogy  of  the  Peierls  resistance  to  a  dislocation 
motion,  is  studied.  They  used  two  models  which  were  subjected  to  a 
"lattice-static"  analysis.  In  the  first  model,  a  one-dimensional 
system,  the  crack  was  depicted  as  two  semi- infinite  chains  consisting 
of  points  of  atoms  linked  horizontally  by  bendable  elements  and  ver¬ 
tically  by  stretchable  elements.  The  chains  were  subjected  to  opening 
forces  applied  vertically  at  the  free  ends .  In  the  two-dimensional 
model  an  infinite  square  lattice  of  atoms  linked  by  stretchable  and 
bendable  elements  was  considered.  The  crack  was  opened  either  by  a 
vertical  wedging  force  or  by  vertical  tensile  forces  distributed  uni¬ 
formly  at  infinity.  The  "lattice-static"  analysis  began  with  an  as¬ 
sumption  concerning  the  form  of  the  atomic  interaction,  and  proceeded 
to  calculate  an  equilibrium  configuration  consistent  with  appropriate 
boundary  conditions.  This  involved  considerable  mathematical  com¬ 
plexity,  and  one  had  to  restrict  the  analysis  to  the  simplest  force 
laws  for  the  linking  elements.  Accordingly  it  was  assumed  that  the 
elements  were  Hookean  up  to  a  critical  breaking  point.  The  lattice 

trapping  behavior  has  been  tested  with  other  atomistic  models,  using 

(9)  (14) 

CMD,  by  P.C.  Gehlen  and  M.F.  Kanninen,  and  A.  Gohar .  These 

results  seemed  to  indicate  that  the  phenomenon  of  lattice  trapping 
is  almost  completely  attenuated  when  more  realistic  interatomic  force 
laws  are  used  in  the  simulation. 

Highly  idealized  modes  have  been  used  in  the  studies  of  Ashurst 
and  Hoover^ 15 ^  and  Weiner  and  Pear.^1^  Ashurst  and  Hoover  used  a 
two-dimensional  triangular  lattice  in  which  the  particles  interact  by 
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truncated  Hooke’s  law  forces  and  the  exterior  stresses  acted  directly 
over  the  free  boundaries.  Their  static  results  for  the  energy,  en¬ 
tropy,  stress  concentration  and  crack  structure  are  consistent  with 
expectations  from  macroscopic  elasticity  theory.  The  dynamic  theory 
of  crack  propagation  is,  according  to  the  authors,  less  well  developed 
and  supersonic  crack  velocities  could  be  produced  by  the  proximity 
of  the  exterior  boundaries.  Under  these  conditions  crack  propagation 
can  outrun  lattice  relaxation.  These  effects  can  be  avoided  if  the 
crack  is  enclosed  in  a  continuum  elastic  medium  in  which  the  stresses 
are  acting  at  infinity.  Then  the  crack  propagation  will  be  only  af¬ 
fected  by  the  natural  relaxation  of  the  strain  field  around  the  crack 
and  will  not  be  influenced  by  the  proximity  of  the  exterior  boundaries. 
This  second  model  is  more  realistic  but  it  introduces  additional  com¬ 
putational  difficulties  that  remain  to  be  solved.  In  the  next  chapter 
we  give  a  full  description  of  the  advantages  and  limitations  of  this 
method,  called  flexible  boundary  conditions. 
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Chapter  Four 


Flexible  Boundary  Conditions 


4.1  Description 

4.1.1  Linear  Elasticity  Solution 

4.1.2  Relaxation  Method 

4.1.3  Green's  Function 

4.2  Computer  Molecular  Dynamics  Determination  of  Elastic  Constants 
in  a  Two-Dimensional  Triangular  Lattice 

4.3  Application  of  Flexible  Boundary  to  a  Two-Dimensional  Crack  Tip 
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4.1  Description 

Several  kinds  of  boundary  conditions  have  been  applied  recently 
to  the  study  of  crack  tip  configurations.  Some  of  them,  periodic 
boundaries,  free  boundaries  with  stresses  applied  directly  on  the  free 
surface,  etc.,  were  summarized  in  Chapter  Three.  In  this  section  we 
give  a  complete  description  of  the  so-called  flexible  boundary  (FB) 
developed  by  Sinclair  et  al.^  The  next  section  describes  the  deter¬ 
mination  of  elastic  constants  for  a  two-dimensional  triangular  lattice, 
a  step  that  is  necessary  in  applying  the  continuum  elasticity  solutions 
at  several  stages  of  FB. 

(2  3) 

Several  methods  of  flexible  boundary  ’  have  been  developed. 

They  differ  only  in  the  way  the  atomic  forces  arising  from  the  inter¬ 
action  of  the  continuum  with  the  atomic  region  are  relaxed.  These 
techniques,  as  noted  by  Sinclair  et  al, ^  have  improved  the  efficiency 
of  the  calculations  by  permitting  the  use  of  smaller  atomic  systems, 
have  yielded  descriptions  of  the  dilatational  effects  arising  from 
non-linearity  in  the  core  of  the  defect,  have  allowed  equilibrium  to 
be  established  between  the  atomic  and  continuum  regions,  and  have  made 

possible  calculations  of  defect  mobility  by  removing  the  constraints 

(2-6) 

imposed  by  rigid  boundaries. 

The  first  step  in  imposing  FB  begins  with  a  perfect  lattice  of 
the  type  of  material  which  is  going  to  be  studied.  The  first  config¬ 
uration  of  the  defect  is  then  generated  with  all  the  atoms  in  the 
system  displaced  according  to  a  linear  continuum  solution.  This  solu¬ 
tion  gives  the  displacement  field  under  fixed  conditions  (type  of  stress 
boundary  conditions,  etc.)  and  the  core  of  the  defect  is  arbitrarily 
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fixed  at  the  center  of  our  system.  The  resulting  array  is  then 
divided  into  three  regions  as  shown  in  Fig.  4-1  for  the  case  of  a 
crack  tip.  Region-I  contains  those  atoms  which  are  going  to  be  moved 
according  to  CMD,  under  the  force  law  of  a  chosen  two-body  central 
potential.  Region-II  is  composed  of  atoms  which  interact  with  atoms 
in  Region-I  or  Region-Ill,  and  contains  all  the  atoms  on  which  a  force 
may  be  exerted  by  at  least  one  atom  in  Region-I.  The  thickness  of 
Region-II  should  thus  equal  the  maximum  range  of  the  interatomic  force 
law.  Region-Ill  is  composed  of  the  exterior  atoms  which  are  always 
displaced  according  to  a  linear  continuum  solution.  The  thickness  of 
Region-Ill  should  correspond  to  the  minimum  distance  to  completely  de¬ 
termine  the  forces  acting  on  each  Region-II  atom  and  is  also  equal  to 
the  maximum  range  of  the  potential  function. 

The  second  step  is  to  relax  Region-I  atoms  using  CMD.  Several 
alternative  procedures  can  be  used  to  reach  a  relaxed  configuration 
in  which  the  forces  on  the  Region-I  atoms  are  less  than  a  predetermined 
value.  The  procedure  followed  here  has  been  explained  in  Chapter  Two; 
it  is  a  method  which  continuously  takes  out  kinetic  energy  from  the 
system  through  a  frictional  force. 

The  third  step  uses  the  residual  forces  acting  on  Region-II 
atoms,  generated  by  the  relaxation  of  Region-I,  to  create  a  displacement 
field  acting  on  the  entire  system.  This  field  is  calculated  by  using 
a  continuum  linear  Green's  function, which  gives  the  displacement 
due  to  a  unit  force  acting  in  a  specified  direction  in  the  presence 
of  the  crack. The  resulting  displacement  field  is  thus  a  sum  of 
Green's  functions  with  origins  at  each  of  the  Region-II  atoms.  As  we 
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Figure  4.1 


Regions  for  the  flexible  boundary 
method  in  the  case  of  a  crack  tip. 
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will  see  in  the  next  sections  the  Green's  function  diverges  logarithmi¬ 
cally  at  the  point  where  the  force  is  applied  so  an  additional  step 
must  be  used  to  move  Region-II  atoms.  In  the  present  procedure  this 
step  is  a  relaxation  of  Region-II  atoms  using  the  frictional  damping 
method  discussed  above. 

If  the  resulting  forces  on  Region-II  atoms  are  not  small  enough, 
the  last  steps  are  iterated  until  all  the  atoms  in  Regions-I  and  II 
have  forces  below  some  chosen  tolerance  level.  A  synopsis  of  the 
total  procedure  constituting  the  flexible  boundary  is  the  following: 

(a)  Generation  of  the  atomic  configuration  of  the  defect 
according  to  a  linear  continuum  solution.  Division 
of  the  system  into  three  regions. 

(b)  Relaxation  of  Region-I  using  CMD. 

(c)  Displacement  of  Regions-I  and  II  according  to  the 
Green's  function  for  the  crack  problem. 

(d)  Relaxation  of  Region-II  atoms  using  CMD. 

(e)  Iteration  of  steps  b,  c,  and  d  until  total  forces 
over  atoms  in  Regions-I  and  II  are  below  a  pre¬ 
scribed  magnitude. 

4.1.1  Crack  Tip  Elastic  Field 

In  proceeding  to  an  analysis  of  the  plane-crack  problem  it  is 
useful  to  distinguish  three  basic  modes  of  crack-surface  displacement. 
Mode  I  corresponds  to  normal  separation  of  the  crack  walls  under  the 
action  of  tensile  stresses;  Mode  II  (sliding  mode)  corresponds  to 
shearing  of  the  crack  walls  in  a  direction  normal  to  the  crack  front; 
Mode  III  (tearing  mode)  corresponds  to  mutual  shearing  parallel  to  the 
crack  front.  Of  the  three  modes,  the  first  one  is  the  most  perti¬ 

nent  to  crack  propagation  in  brittle  solids  in  which  we  can  visualize 
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the  crack  extension  by  progressive  stretching  and  rupture  of  cohesive 
bonds’  across  the  crack-plane. 

In  our  study  we  will  concentrate  on  Mode  I  fracture  where  small 
plasticity,  with  possibility  of  formation  of  dislocations,  can  result 
around  the  crack  tip.  At  greater  distances  from  the  crack  tip  the 
linear  elasticity  theory  can  be  applied.  For  non-brittle  solids, 
blunted  cracks  with  high  production  of  dislocations  and  large  plastic 
regions  around  the  crack  tip  or  completely  surrounding  the  crack  can 
be  produced. 

In  this  section  we  give  the  final  solution  for  the  linear  elas¬ 
ticity  crack  tip  field,  stresses  and  displacements  in  the  Mode  I 
fracture,  which  are  used  in  the  first  step  of  the  FB  technique.  The 
standard  solution  to  this  problem  involves  searching  for  a  suitable 
"stress  function"  which  satisfies  the  "biharmonic  equation"  of  linear 
elasticity  theory  (fourth-order  differential  equation  including  the 
condition  of  equilibrium, compatibility  of  strains,  and  Hooke’s  law), 
in  accordance  with  appropriate  boundary  conditions.  ^0*11)  The  com_ 
ponents  of  stress  and  displacements  are  then  determined  directly  from 
the  stress  function. 

An  analytical  technique  has  been  developed  by  Westergaard, 

Muskhelishvili  and  others  for  the  special  case  of  plane-crack 
(12  13) 

geometry.  ’  The  model  is  the  "sharp  crack"  approximation,  where 
it  is  assumed  that  the  crack  tip  in  the  unstressed  state  is  perfectly 
sharp  and  that  the  crack  walls  remain  free  of  stresses. 

The  solutions  for  the  field  near  the  crack  tip  (i.e.,  at  dis¬ 
tances  from  the  crack  tip  small  compared  to  the  length  of  the  crack) 
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take  the  following  analytic  form  for  the  isotropic  case 
Mode  I 


(10) 


(Fig.  4.2): 
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where  K  =  (3->0/(l+>0  for  plane  stress,  K  =  (3-4  >0  for  plane  strain, 

E  is  the  Young’s  modulus  and  V>  the  Poisson's  ratio.  The  quantity 
is  called  the  stress-intensity  factor;  it  is  dependent  on  the  boundary 
conditions  of  the  crack  system.  In  this  case  its  value  is: 

=  rry^^  0*^ 

=  tensile  stress.  Mode  I,  applied  at  infinity 
Z  c  =  crack  length 

We  can  reduce  these  expressions  to  the  simple  form: 
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The  stress  intensity  factor  depends  only  on  the  applied  stress 
and  the  crack  geometry;  it  determines  the  intensity  of  the  local  field. 
The  other  factors  depend  on  the  spatial  coordinates  about  the  tip,  and 
determine  the  distribution  of  the  field. 

4.1.2  Atomic  Relaxation 

The  next  step  of  the  flexible  boundary  technique  consists  in 
the  relaxation  of  Region-I  toward  a  configuration  of  very  low  forces 
with  Region-II  and  Region-Ill  atoms  held  fixed.  As  indicated  before, 
the  procedure  followed  here  involves  a  frictional  damping  applied  to 
the  "equations  of  motion"  after  a  certain  period  of  normal  simulation. 
The  final  temperature,  after  application  of  frictional  damping,  is 
reduced  to  approximately  3°K.  Further  damping  to  still  lower  tempera¬ 
tures  does  not  lead  to  any  significantly  different  configuration  of 
the  relaxed  atomic  positions. 

4.1.3  Green* s  Function 

The  next  step  of  our  procedure  which  is  different  from  the 
previous  flexible  boundary  methods  is  to  generate  a  new  displacement 
field  based  on  the  residual  forces  on  Region-II  atoms.  These  forces 
arise  from  the  constraint  imposed  by  Region-Ill,  composed  of  atoms 
displaced  according  to  a  continuum  linear  elasticity  solution,  over 
Region-I . 

The  displacement  field  necessary  to  cancel  these  forces  is  calcu¬ 
lated  by  a  sum  of  Green's  functions  corresponding  to  each  force  on 
Region-II  atoms.  The  Green's  function^  gives  the  displacement  field 
caused  by  a  unit  point  force  applied  in  a  certain  direction  at  some 
distance  from  the  crack.  This  Green's  function  must  include  the 
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influence  of  the  crack  surfaces  on  which  the  boundary  conditions 

must  be  satisfied  (no  forces  can  act  normal  to  a  free  surface). 

The  resulting  displacement  field  due  to  the  Green's  function 

(8) 

has  been  calculated  by  Hirth  et  al.  for  the  case  of  isotropic 
materials.  The  displacement  field  due  to  the  image  forces  which 
cancel  the  normal  stresses  on  the  two  crack  surfaces  is: 
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where  and  are  the  point  forces  in  the  x  and  y  directions, 

=  F^/F2»  K  =  3-4 v  for  plane  strain,  y.  is  the  shear  modulus,  V  ,  the 
Poisson's  ration  and 

Or  1+  Pl  +  IP  wl  9  * 

2. 

£  =  I  +  P1upw  ti± 

2. 

f  ?  on  c&/tj 
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r,  ©  and  p ,  'f  are  the  cylindrical  coordinates  defining  the  positions 
of  the  force  source  and  point  where  the  displacement  is  evaluated. 

The  field  of  the  forces  themselves  that  gives  the  total  field 
when  added  to  the  displacement  field  of  the  image  forces  has  been 
derived  by  Hirth^  and  Eshelly.  The  displacement  field  is: 
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We  rederived  these  results  following  the  procedure  indicated  by 
Hirth. ^  An  elastic  and  isotropic  Green’s  function  is  calculated  in 
the  same  way  as  for  the  evaluation  of  potential  fields  in  electro¬ 
statics.  In  this  case,  there  is  a  correspondence  between  the  elastic 
displacement  <x  and  the  electrostatic  potential  V  and  between  the 
sources  generating  the  fields,  force  f,  and  the  charge  density  p .  One 
finds  that  the  ith  component  of  the  displacement  iA.  ■  ■  [  f-  r‘)  caused  by  a 
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unit  point  force  <^(r')  applied  in  the  jth  direction  at  the  point  rT 
(7) 


is 


J-  {$.,  j  )  )  (4-5) 

\vjl\  A+Z^\  t>  J  J  K  y 


(T-ir^)  is  called  the  Green's  function  for  the  elastic  displacements 
of  an  isotropic  system  withoLline.  constants  X  and/i.  A  continuous 
distribution  of  forces  in  an  elastic  medium  causes  displacements 

y 

Equation  (4-5)  gives  the  response  of  an  infinite  body  to  a  point 
force.  In  a  finite  body,  boundary  conditions  must  be  satisfied.  The 
displacements  in  a  finite  body  subjected  to  a  point  force  can  be  de¬ 
scribed  as  a  superposition  of  the  displacements  (4-5)  and  displace¬ 
ments  caused  by  "image"  forces  applied  on  the  external  surface  of  the 
body  in  order  to  satisfy  boundary  conditions. 

In  our  case  we  have  rederived  Eqs.  (4— 4b)  applying  (4-5)  and 
(4-6)  to  the  case  of  two  line  forces  %  kx’J  Scy'J  and  fy  liy‘)  . 

The  forces  are  uniformly  distributed  along  a  line  parallel  to  the  Z 

axis.  ¥  and  ¥  are  the  forces  per  unit  length  directed  in  the  x  and 
x  y 

y  direction  respectively. 

The  method  we  have  applied  is  a  powerful  procedure  for  the  solu¬ 
tion  of  problems  in  the  continuum  theory  of  elasticity.  For  example, 
the  displacement  field  due  to  point  and  line  defects  (interstitials, 
dislocations,  etc.)  can  be  easily  determined  if  the  response  of  the 


54 


body  to  a  point  force  (4-4)  or  a  line  force  (4-4b)  is  known.  The 
interstitial  solution  would  be  given  by  the  displacement  field  due 
to  three  perpendicular  double  forces,  and  an  edge  dislocation  by  the 
displacement  field  of  two  line  forces  coupled  without  moment. 

To  make  an  estimate  of  the  range  of  validity  of  the  continuum 
Green’s  function,  a  comparative  study  has  been  done  between  the  elastic 
displacement  provided  by  this  function  and  the  displacement  obtained 
using  CMD  at  short  distances  from  the  point  force.  The  CMD  experiment 
consisted  of  applying  an  external  unit  point  force  at  the  center  of  a 
perfect  lattice  and  calculating  the  final  displacement  field.  The 
external  force  was  applied  to  an  atom  surrounded  by  an  infinite  lattice 
with  the  boundary  atoms  located  according  to  the  continuum  Green's 
function  (Eq.  (4-46).  The  final  CMD  displacement  field  gave  values 
about  2%  less  than  the  Green's  function  for  the  first  nearest  neighbors 
while  at  greater  distances  from  the  source  the  discrepancy  was  much 
less.  Therefore,  the  lattice  Green's  function  can  be  applied  at  dis¬ 
tances  greater  than  or  equal  to  the  nearest  neighbor  separation.  The 
asymmetry  of  the  displacement  field  in  the  CMD  case  is  due  to  the  non¬ 
linear  effects  of  the  potential  function. 
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Table  4.1  Atomic  displacements  (A)  given  by  the  Green’s 

function  of  a  unit  point  force  directed  in  the 
positive  x  direction  and  CMD  simulation.  The 
displacements  have  been  measured  at  several 
interatomic  spacings  from  the  force  along  the 
x  axis.  The  interatomic  spacing  is  4.1290A. 
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4.2  Computer  Molecular  Dynamics  Determination  of  Elastic  Constants 
in  a  Two-Dimensional  Triangular  Lattice 

The  macroscopic  elastic  constants  of  a  2-D  triangular  lattice 
have  been  determined  by  CMD.  The  system  is  composed  of  rare  gas  atoms 
interacting  with  nearest  neighbors  through  a  Lennard— Jones  potential. 
The  CMD  simulation  was  conducted  by  applying  different  kinds  of  stress 
on  the  free  surfaces  of  the  system,  which  is  composed  of  an  8x8  array 


of  particles. 

The  elastic  constants  were  determined  in  the  limit  of  small 
deformations  where  the  linear  elasticity  theory  can  be  applied.  A 
frictional  damping  procedure  was  used  during  the  atomic  relaxation. 

The  final  position  of  static  equilibrium  corresponds  to  a  state  of 
minimum  total  potential  energy  and  its  variation  during  the  experiment 
must  be  equal  to  the  work  done  by  the  exterior  forces. 

The  present  model  is  similar  to  the  one  studied  by  Ashurst  and 
Hoover . They  have  assumed  a  two-dimensional  triangular  crystal 
in  which  particles  interact  with  nearest  neighbors  springs  of  elastic 
constant  K.  Their  calculations  showed  the  isotropic  form  of  the 
stress-strain  relations  and  thec^On^e  constants  were:  ^  ^  /  4 

(See  Eq.  (4-11).) 

The  validity  of  the  Ashurst-Hoover  analytical  results  was  tested 
independently  by  our  CMD  experiments.  Both  cases  gave  the  same  results 
for  Is:  f  ,  where  the  equilibrium  position  corresponds 

to  the  interatomic  distance  between  nearest  particles.  Logically,  for 


small  displacements,  any  potential  would  give  the  same  elastic  con¬ 
's4-  \j  i 

stants  if  its  curvature  at  the  position  of  equilibrium  <x  _ — - 

O  V4  I  wn 


is  fitted  to  the  same  value. 
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By  linear  elasticity  theory, ^18, 19^  the  elastic  constants  in  a 
two-dimensional  crystal  are  completely  specified  if  the  tensor  T  is 
known: 

r=  t  & 

in  two  dimensions 
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where  5  cj*  <r*y  '  °"V* 

and  c*‘  =  Cjl-  =C|>  r  Ci*  =  ° 

The  calculations  of  the  four  constants,  C^,  C 22  anc*  ^33’  were 

performed  in  four  separate  CMD  simulations. 

(A)  Simulation  1  -  Calculation  of  C22‘ 

In  this  calculation  we  began  with  the  equation 

:  ci.  +•  L  yy 

We  made  €xx:o  by  ensuring  no  atomic  displacements  in  the  x  direction, 
then 

5~yy  -  c*2.  £  yy 

where  ry. :  [H ,  <=yy :  *■'*•«/!-„  ;f  is  the  exterior  force  acting  on  each 
atom,  d,  the  nearest  neighbor  interatomic  spacing  and  L,Lothe  final 
and  original  length  of  the  system  respectively.  C^2  is  determined 
from  the  slope  of  the  curve  versus  €yya t  the  origin.  Strains 
were  also  calculated  by  measuring  the  relative  deformations  within 
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the  system.  In  this  way  the  boundary  effects  were  avoided. 

(B)  Simulation  2  -  Determination  of 

<Tyy  =  Cl,  eXJ«  f  €yy 


Now,  we  let  contractions  of  the  system  in  the  x-direction  while  apply¬ 
ing  a  stress  <Tyy .  The  two  strains  <E and  £yy  were  measured  and  the 
value  of  C^2  was  known  from  Simulation  1. 

(C)  The  same  procedure  was  employed  to  calculate  C^,  but  now 
applying  a  stress  in  the  x-direction.  The  symmetry  of  the  T  tensor 
was  checked  by  determining  C^-  As  expected  from  the  isotropic  form 
of  T  (Eq.  (4-11)),  the  values  found  for  and  C ^  resulted  equal  to 
C22  and  C22  respectively. 

(D)  Finally,  to  determine  a  shear  stress, 

tangent  force  to  the  free  surface,  was  applied  on  each  surface  atom, 
whereby  definition,  £ xy  =  8^  9  .  The  resulting  elastic  coefficients 

for  our  2-D  lattice  are 


C„  z  Cu  :  l  -IS 6  .  /6  ity'wei  / Om 
C(A  7  Clt  :  418.72.  /C/m 

C :  $3 7.  <*4  Ayvinlovn 


If  we  express  Eq.  (4-7)  in  the  form: 
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the  value  of  the  o^cuwe  constants  X  and  }x  can  be  determined  using  (4-10) 
X~  }Jl  -  4IX.72.  dy'nfc*  /  c-wi 


We  can  see  that  these  values,  determined  by  CMD,  coincide  with  the 
analytical  calculations  of  Ashurst  and  Hoover,  assigning  to  the  con¬ 
stant  K  the  value  fLMbl 1_  [  =  <?£7 

5  *  (^luUtWHwl 

The  determination  of  the  elastic  constants  E,  Young’s  modulus, 
andy,  Poisson’s  ratio,  was  carried  out  by  inverting  the  tensor  T. 
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In  the  following  CMD  experiments  we  have  used  the  elastic  con¬ 
stants  derived  in  this  section.  When  it  was  necessary  to  use  elastic 
continuum  solutions  in  three  dimensions,  for  example  in  step-1  of  FB/ 
where  the  initial  configuration  of  the  defect  was  generated,  and  in 
step-3  where  Green's  function  was  applied,  we  have  used  the  equivalence 

between  our  2-D  case  and  the  3-D  plane  strain  with  the  same  <AdWne 
(16) 


constants. 
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Surface  energy  is  another  physical  constant  of  great  importance 
in  mechanics  of  fracture.  In  our  case  we  need  it  to  determine  the 
theoretical  Griffith’s  critical  stress.  The  surface  energy  is  taken 
as  the  work  required  to  completely  separate  the  crystal  along  some 
given  plane  starting  from  the  undeformed  perfect  configuration. 

When  the  atomic  forces  are  specified  by  a  potential  function  and 
there  is  a  finite  cut-off  distance,  the  surface  energy  is  given  by  the 
difference  of  energy  on  the  potential  curve  between  the  static  position 
of  equilibrium  and  the  cut-off  distance  which  would  correspond  to  the 
energy  to  break  a  bond  times  the  number  of  bonds  per  unit  area  along 
the  plane  of  crack  propagation. 

We  have  determined  by  using  this  procedure  the  surface  energy  of 
the  Lennard-Jones  system  and  three  Morse  potentials.  The  values  are 
listed  in  Table  5.1  of  the  next  chapter. 

4.3  Application  of  Flexible  Boundary  to  a  Two-Dimensional  Crack  Tip 

The  procedure  described  in  Section  4.1  has  been  applied  to  a  two- 
dimensional  triangular  lattice  whose  atoms  interact  through  a  Lennard- 
Jones  potential.  Fig.  4.3. 

To  demonstrate  the  efficiency  of  the  flexible  boundary  technique 
as  applied  to  small  systems,  several  tests  were  carried  out  with  dif¬ 
ferent  numbers  of  particles  and  two  different  stresses.  Figure  4.4a 
shows  the  variation  of  the  crack  tip  bond  length  with  the  size  of  the 
system  for  a  stress  intensity  factor  close  to  the  Griffith's  critical 
value,  which  has  been  determined  using  the  elastic  constants  and  sur¬ 
face  energy  of  the  system  (Section  4.2) 

(r  •  VrrtJ=  VIeJ  =  •  i°’z  Jy*15  '  ^ 


Figure  4.4  a,  -  &  Variation  of  crack  tip  bond  length  with  system  size,  [axb]  denotes  the 

number  of  particles  along  x  and  y  direction  in  Region-I.  Results  for 
fixed  boundary  condition  (0),  flexible  boundary  condition  (•) .  (a) 

cr^Tc=3.17‘l0-2  dynes  cm~i/2.  (b)  <r  \/^ c=3 . 17 *  10-^  dynes  cm-^'^. 
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The  high  efficiency  of  the  procedure  is  shown  in  Fig.  4.4b  where 
a  small  stress  intensity  factor  was  applied  to  ensure  that  the  linear 
elasticity  theory  is  valid  at  any  point  of  the  system.  The  configura¬ 
tion  of  equilibrium  using  FB  could  be  sufficiently  accurately  determined 
with  a  small  system  composed  of  only  84  particles  whereas  considerably 
larger  systems  were  required  to  obtain  the  same  crack  tip  bond  length 
using  fixed  boundary  condition. 

At  high  stresses,  where  a  non-linear  region  exists  around  the 
crack  tip,  the  flexible  boundary  procedure  is  not  any  more  efficient, 
as  Fig.  4.4a  shows  the  results  given  by  FB  vary  with  the  size  of  the 
system  about  as  much  as  the  results  obtained  with  fixed  boundary.  The 
reason  for  this  behavior  is  probably  that  the  elastic  continuum 
Green's  function,  which  is  responsible  for  the  final  displacement  of 
the  boundaries,  has  been  derived  assuming  the  system  is  entirely  com¬ 
posed  of  a  linear  elastic  material  and  does  not  consider  the  non-linear 
properties  of  the  region  surrounding  the  crack  tip.  This  effect  is 
clearly  more  important  at  high  stresses  and  when  the  system  size  is  no 
longer  large  compared  to  the  size  of  the  non-linear  region.  In  this 
case  it  will  be  necessary  to  work  with  larger  systems  in  order  to 
obtain  the  correct  crack  tip  configuration. 

Table  4.2  shows  the  variations  with  system  size  of  the  maximum 
force  acting  on  atoms  in  Region-I  and  Region-II.  Also  given  are  the 
crack  tip  bond  lengths  obtained  with  fixed  and  flexible  boundaries. 

If  the  forces  acting  on  Region-II  atoms  are  zero  after  the  re¬ 
laxation  of  Region-I,  we  can  assume  that  the  continuum  linear  elasticity 
theory  is  valid  in  this  region  and  that  the  boundaries  do  not  impose 
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any  constraint  on  the  crack  tip.  Unfortunately,  the  magnitude  of 
these  forces  decreases  very  slowly  by  increasing  the  number  of  particles 
(fourth  column  of  Table  4-2)  being  necessary  to  use  very  large  systems 
to  assume  that  the  continuum  solution  at  the  boundaries  is  valid. 
Flexible  boundaries  provide  a  procedure  to  cancel  these  forces  for 
small  systems  and  as  we  have  seen  its  efficiency  only  decreases  in 
the  presence  of  large  non-linear  regions  around  the  crack  tip. 

Figures  4-5  to  4-16  show  the  strain  fields  and  rotations  around 
the  crack  tip  at  different  stages  of  the  computations:  (a)  initial 
step  or  continuum  elasticity  solution,  (b)  after  relaxation  of  Region-I, 
fixed  boundary  and  (c)  after  the  application  of  FB.  The  principal 
strains  and  rotations  have  been  calculated  using  the  procedure  given 
in  the  next  chapter. 

The  difference  between  fixed  and  flexible  boundaries  is  shown 

in  Figs.  4.5a,  b  and  c,  where  a  stress  intensity  factor  trVrrc  equal  to 

_ vfchto  _ .♦**»« 

1.13(0*VtrcV  was  applied.  The  critical  value  (<rvn corresponds 

to  the  theoretical  value  derived  by  Griffith  and  was  calculated  from 

the  elastic  constants  and  surface  energy  of  the  system  . 

Figure  4.5b  shows  the  crack  tip  configuration  after  the  relaxa¬ 
tion  of  Region-I,  fixed  boundary.  The  crack  tip  bond  length  has  not 
reached  the  maximal  permissible  length,  which  corresponds  to  the 
cut-off  of  the  potential,  so  one  may  conclude  on  this  basis  that  the 
critical  value  of  the  stress  intensity  factor  is  greater  than 

; _ ,H>eo 

1.  lSCQvnc)  .  However,  as  shown  in  Fig.  4.5  c,  which  corresponds  to 
a  configuration  obtained  with  FB,  the  crack  tip  bond  has  already 
reached  the  maximum  bond  length  at  a  stress  intensity  factor  of 
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l-lSC^l^nc)^  *fl,  and  the  crack  is  propagating  by  breaking  bonds. 
Therefore,  the  critical  stress  intensity  factor  is  actually  less  than 
i.l3(rl/ic  )*“ 

Figures  4.6  and  4.7  show  the  principal  strain  fields  and 
rotations  around  the  crack  tip  for  the  different  boundaries.  The 
linear  elasticity  solution,  Eq.  (4- la),  yields  a  state  of  biaxial 
tensile  strain  in  front  of  the  crack  tip  and  higher  values  of  shear 
deformation  at  60°  and  120°  from  the  positive  x  axis  (see  Fig.  4.6a). 
The  non-linear  effects  given  by  CMD  simulations  are  shown  in  Figs. 

4.6b  and  4.6c,  which  correspond  to  the  amplifications  of  these  fields 
in  the  immediate  vicinity  of  the  crack  tip.  The  brittleness  of  the 
material  has  not  changed  substantially  the  distribution  of  the  strain 
around  the  crack  tip  and  the  only  remarkable  characteristics  are  the 
higher  values  of  the  strains  when  the  FB  procedure  was  applied. 

The  absence  of  dislocation  in  this  material  which  would  indicate 
a  brittle  behavior  can  be  explained  by  applying  the  Kelly,  Tyson  and 
Cottrell  criterion  (KTC)  discussed  in  Chapter  Three.  This  criterion 
for  brittle  behavior  is 
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where  the  second  term  depends  on  the  maximum  values  attained  at  the 
crack  tip.  Its  value  according  to  linear  elasticity  theory  is  0.5 
for  mode  -  1  fracture. 


The  ideal  shear  and  cohesive  stresses  were  determined  by  two 
separate  simulations,  where  a  shear  and  a  biaxial  tensile  stress  were 
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applied  to  a  perfect  lattice  (Figs.  4.10  and  4.11).  We  will  show  in 

the  next  chapter  that  the  ideal  shear  stress  determined  in  this  way  is 

(21) 

an  overestimate.  As  proposed  by  Orowan  we  have  taken  for  the 
ideal  shear  stress  the  value  determined  in  a  simulation  where  the 
atomic  sliding  in  the  x  direction  is  accompanied  by  a  small  displace¬ 
ment  in  the  positive  y-direction  (Fig.  3.11,  dashed  line).  The  rela¬ 
tion  between  the  ideal  shear  and  cohesive  stress  is  0.83,  implying 
that  the  system  should  show  brittle  behavior. 

Finally,  it  has  been  found  that  the  stress  intensity  factor 
necessary  to  break  the  crack  tip  bonds  and  propagate  the  crack  is 
higher  than  predicted  by  Griffith’s  theory. 

C  M  D  IS  •  I  O'*  e/y  ^  Cl-o) 

Gri  f  f C  tVj  Hieovy  (<r“  1/ttc  J  c  5  !0  1  •  c**i  <•  (.2-0 J 

The  CMD  critical  stress  intensity  factor  was  determined  by  ap¬ 
plying  successively  higher  values  to  <rVr»c.  until  crack  propagation  was 
obtained. 

The  critical  value  depends  on  the  atomic  model  used  in  the 

simulation.  In  our  case  we  have  assumed  no  interaction  between  atoms 

situated  on  different  crack  surfaces.  If  interaction  is  possible  when 

the  distance  between  a  pair  of  atoms  is  less  than  the  cut-off  range 

of  the  potential  function  the  stress  intensity  factor  can  be  higher 

than  the  value  obtained  above.  Similar  results  were  found  by  Gehlen 
(22) 

et  al.  in  the  case  of d -iron,  the  critical  stress  Intensity  factor 
being  about  three  times  the  Griffith  value.  These  results  are  not 
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surprising  for  several  reasons.  First,  Griffith's  criterion  is  con¬ 
sidered  a  necessary  condition  but  not  sufficient  for  crack  propagation. 
Secondly,  Griffith's  equations  are  based  upon  macroscopic  linear 
elasticity  continuum  properties.  They  do  not  consider  energy  balance 
on  a  local,  atomic  scale  or  generation  of  dislocation  at  the  crack 
tip,  although  the  last  possibility  was  not  observed  in  this  material. 
Further  discussion  of  these  results  and  some  of  the  recently  developed 
atomic  models  which  could  explain  this  behavior  is  given  in  the  last 
chapter. 
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Figure  4.5a  Initial  atomic  configuration  based  on  linear 

elasticity  theory  of  a  two-dimensional  crack 
embedded  in  an  infinite  medium.  Sytem  consists 
of  436  particles  arranged  in  a  triangular  lat¬ 
tice.  The  crack  tip  bond  is  shown  by  the  line 
segment  A,  its  length  is  less  than  the  cut-off 
of  the  potential  function. 
ry"rFc=1.13(  'F  l/irc)£heoretical. 
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Same  system  as  Fig.  4.5a  after  the  relaxation 
of  Region-I  (fixed  boundary) .  The  crack  tip 
bond  is  A1.  Its  length  is  now  10%  greater  than 
A  (Fig.  4.5a)  but  still  less  than  the  cut-off 
range  of  the  potential.  However,  relaxation 
of  the  forces  acting  on  the  particles  at  the 
boundaries,  by  using  FB,  reveals  that  the 
critical  stress  intensity  factor  has  been 
exceeded  (see  next  figure) . 


Figure  4.5b 
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Figure  4.5c  Atomic  configuration  using  flexible  boundary. 

The  crack  starts  to  propagate  by  breaking 
bonds  because  the  critical  stress  intensity 
factor  was  exceeded.  The  crack  tip  bond  is 
now  A".  Previous  crack  tip  bonds  have  been 
broken  (dashed  lines). 
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Principal  strains  at  the  initial  configuration 
(linear  elasticity  theory)  J  dilatation  strain, 
|  compressive  strain. 


Figure  4.6a 
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Figure  4.6b  Principal  strains  after  the  first  relaxation 

(fixed  boundary)  $  dilatation  strain, 

|  compressive  strain. 
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Figure  4.6c  Principal  strains  in  the  final  relaxed 

configuration  (flexible  boundary). 

$  dilatation  strain,  |  compressive 
strain. 
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Figure  4.7a  Rotations  at  the  initial  configuration 

(linear  elasticity  theory)  f  indicates 
magnitude  and  direction  of  the  rotation. 
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Figure  4.7b 


Rotations  after  the  first  relaxation 
(fixed  boundary)  f  indicates  magnitude 
and  direction  of  the  rotation.  ■ 
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Rotations  after  applying  flexible  boundary,  f 
indicates  magnitude  and  direction  of  the 
rotation. 


Figure  4.7c 
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Figure  4.8a  Principal  strains  at  the  initial  configuration 

(linear  elasticity  theory)  J  dilatation  strain, 
|  compressive  strain. 
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Figure  4.8b  Principal  strains  after  the  initial  relaxation 

(fixed  boundary)  £  dilatation  strain, 

J  compressive  strain. 
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Figure  4.8c 


Principal  strains  (flexible  boundary) 
dilatation  strain  j  compressive  strain. 


81 


i  <1  i  1  i  *  i 

A  A  1  1  1  l  l 

WiV.v. 

'/t'tVt''  ■  • 

T  t  T  T  t  t  t  . 

T  T  T  t  t  t  t  r 

t  r  t  t  t  t  t 


i 


Rotation  field  at  the  initial  configuration 
(linear  elasticity  theory),  f  indicates 
magnitude  and  direction  of  rotation. 


Figure  4.9a 
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Figure  4. 9b  Rotation  field.  Fixed  boundary. 

f  indicates  magnitude  and  direction 
of  rotation. 
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Figure  4.9c 


Rotation  field.  Flexible  boundary. 

f  indicates  magnitude  and  direction 
of  rotation. 


Biaxial  stress  versus  strain  in  the  case  of 
uniform  expansion  of  a  perfect  lattice  whose 
atoms  interact  through  a  Lennard-Jones  poten 
tial  (Kr) .  Point  a  on  the  curve  corresponds 
to  the  ideal  cohesive  stress,  point  b  is  the 
rupture  point,  where  all  the  atomic  bonds  re 
the  cut-off  distance  of  the  potential  functi 


85 


TRIANGULAR  LATTICE 


DISPLACETOTTy'ATOrtlC  SPACENENT 

Figure  4.11  Shear  stress  versus  displacement  for  the  shearing  of 

two  rows  of  atoms  which  interact  through  a  Lennard- 
Jones  potential  (Kr) .  Curve  a  was  obtained  with  dis¬ 
placements  along  the  x  direction.  Curve  b  permitted 
simultaneous  relaxation  along  the  y  direction. 
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Chapter  Five 

Computer  Molecular  Dynamics  Studies  of 
Plastic  Deformation  around  a  Crack  Tip 


5 . 1  Introduction 

5.2  Ductile  and  Brittle  Behavior 

5.3  Computer  Molecular  Dynamics  Studies  of  Ductile  versus 
Brittle  Behavior 
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5 . 1  Introduction 

Two  continuum  treatements  to  predict  brittle  or  ductile  behavior 
from  known  material  properties  have  been  studied  using  CMD.  In  ad¬ 
dition  to  these  two  tests,  the  strain  fields  due  to  several  atomistic 
processes  of  plastic  relaxation  around  the  crack  tip  have  been 
determined. 

The  first  section  of  this  chapter  summarizes  the  Rice  and 
Thomson  (RT) ^  and  Kelly,  Tyson  and  Cottrell  (KTC) ^  criteria  for 
predicting  brittle  or  ductile  behavior.  Section  two  describes  the 
atomistic  model  used  in  the  simulation.  Several  properties  of  the 
material,  elastic  constants  and  ideal  stresses,  were  determined  in 
additional  tests  by  applying  CMD  to  a  block  of  perfect  lattice  under 
different  states  of  stress.  The  continuum  criteria  to  predict  brittle 
or  ductile  behavior,  KTC  and  RT,were  studied  and  compared  with  atomis¬ 
tic  preductions  given  by  CMD.  Finally  the  rotations  and  strain  field 
due  to  plastic  deformation  around  the  crack  tip  were  studied  and 
analyzed  numerically. 

5 . 2  Ductile  and  Brittle  Behavior 


Two  criteria  have  appeared  recently  to  predict  brittle  or  ductile 

behavior  in  a  known  material  under  certain  states  of  stress.  Kelly, 

Tyson  and  Cottrell  proposed  that  brittle  fracture  would  be  observed 

if  / <r )  ,  ,  >  ,  where  the  subscripts  "ideal"  and  "max" 

v  c' ideal  c  max’ 

refer  respectively  to  the  ideal  properties  of  a  perfect  lattice  and 
to  the  maximum  values  attained  at  a  crack  tip.  is  calcu- 

LUcAA 

lated  using  linear  elasticity  theory  (for  mode-1  fracture  it  is  ~  0.5) 

(1) 


and  «n  /!r.)ldeal 


is  a  property  of  the  material.  Rice  and  Thomson 
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on  the  other  hand,  assumed  that  a  necessary  criterion  for  brittle 
fracture  is  stability  against  the  emission  of  dislocations  from  the 
crack  tip  and  have  derived  an  analytical  expression  using  linear 
elasticity  and  the  Peierls  model  of  a  dislocation  core.  They  used 
two  different  treatments  to  predict  the  type  of  fracture.  The  first 
one  is  based  on  a  three-dimensional  model  which  calculates  the  energy 
to  activate  a  dislocation  loop.  The  second  one  uses  a  two-dimensional 
model  with  the  dislocation  emitted  from  the  crack  tip  lying  parallel 
to  the  crack  front.  In  both  cases,  three  forces  are  supposed  to  be 
acting  on  the  dislocation  near  the  crack  tip:  (a)  the  force  due  to 
the  stress  field  surrounding  the  crack,  (b)  the  surface  tension  force 
caused  by  creating  more  surface  at  the  blunted  crack,  and  (c)  the 
image  force  of  the  dislocation  in  the  free  surface  of  the  crack.  The 
first  force  repels  the  dislocation,  and  the  other  two  attract  it  toward 
the  crack  tip,  generating  a  position  of  unstable  equilibrium  where  the 
resultant  force  is  zero.  Both  treatments,  nevertheless,  produced 
similar  quantitative  results and  we  have  studied  the  treatment 
based  on  the  two-dimensional  model  because  it  is  more  closely  related 
to  our  simulation  setup. 

Figure  5.1  shows  the  geometry  of  the  dislocation  and  the  crack 
tip.  The  dislocation  has  a  Burgers  vector  perpendicular  to  the 
dislocation  line  (edge  dislocation);  it  is  emitted  at  an  angle  of  60° 
from  the  positive  x  axis  (direction  of  maximum  shear  stress  in  mode-1 
fracture).  Under  this  condition  the  force  acting  on  the  dislocation 
due  to  the  stress  field  produced  by  the  crack  tip  in  mode-1  fracture 
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Figure  5.1  Dislocation  and  crack  configuration.  The 

Burgers  vector  ,  perpendicular  to  the 
dislocation  line,  corresponds  to  an  edge 
dislocation  emitted  by  the  crack  tip.  Rice- 
Thomson  criterion  assumes  that  the  crack  tip 
will  be  blunted  by  emission  of  dislocations 
(plastic  behavior)  when  the  core  of  the  dis¬ 
location  §0is  greater  than  a  critical  dis¬ 
tance,  Jjt.  This  critical  distance  corresponds 
to  a  position  of  unstable  equilibrium.  For 
Sj  <  ^cthe  dislocation  is  attracted  by  the 
crack  tip  and  for  £  >  ^  it  is  repelled. 


9Q 


I 


{jrcty-  L  J 


where  ^  :  f/|i  ,  <)>:  60° y  y^  -  Ji/m  j>  Ctrl,  [^H-)  5IH , 

c  are  the 

exterior  stress  and  half  crack  length  respectively. 

(3) 

The  attractive  image  force  is 

L.  ,  ,  f  (5-*) 

<5 


where  E  is  the  Young's  modulus  and  )>  the  Poisson's  ratio. 

The  third  force  is  produced  by  the  surface  tension  when  the 
crack  is  blunted  (increment  of  free  surface)  and  is  given  by 


Ot 


■nr 


{s-i) 


Vi 


S.n  , 


^  is  the  surface  energy,  <*  r  JL 

The  radius  of  the  dislocation  core,  Sjd  ,  has  been  determined  using 

(3) 

the  Peierls  model. 

The  attractive  forces,  ^  and  fg,  proportional  to  £  '  and 
respectively,  are  greater  than  when  ^  is  small,  and  the  opposite 
occurs  when  {j  is  large.  Hence,  if  the  dislocation  is  able  to  sur¬ 
mount  the  position  of  unstable  equilibrium,  ,  it  will  be  driven 
away  until  it  reaches  some  obstacle  or  it  cannot  overcome  the  resis¬ 


tance  of  the  lattice. 
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The  critical  distance,  ,  at  which  the  dislocation  reaches  an 
unstable  equilibrium  is,  f rom  Eqs.  (5-1),  (5-2)  and  (5-3) 

f  _  g-  v/in if  ^  V4J_  -  t-k  J _  -  'tAJC  —  =  o  C5_ZyJ 

>  total  ~  Vfci x\cJ  H  <iTT  £fc(j-  Tr 

after  replacing  E  by  2)x(l+V)  and  making  ^  =  6  0°. 

At  this  point,  the  Rice-Thomson  criterion  assumes  that  if  the 
position  of  equilibrium  ^c,  obtained  from  Eq.  (5—4),  is  less  than  the 
dislocation  core,  £  ,  spontaneous  generation  is  possible  and  the  crack 
suffers  ductile  fracture. 

Rice  and  Thomson  have  determined  the  dislocation  core  and 
the  critical  distance  £  for  a  wide  range  of  materials,  and  based  on 
these  results  they  found  that  a  sharp  cleavage  crack  is  stable  in  a 
wide  range  of  crystal  types.  In  general,  face-centered  cubic  materials 
are  able  to  emit  dislocations.  Ionic  and  covalent  crystals  are  stable 
against  dislocation  emission,  and  the  body-centered  cubic  crystals  are 
intermediate  between  brittle  and  ductile  materials. 

5.3  Computer  Molecular  Dynamics  Studies  of  Ductile  versus  Brittle 

Behavior 

The  simulation  system  we  have  studied  is  composed  of  436  part¬ 
icles  forming  a  triangular  lattice  and  interacting  through  a  Morse 
potential  (Fig.  5.2) 

V(V,)r  ol  Voj)  -  i  'r°-0  j 


1 


1.4 


1.8 


Figure  5 . 2 


r 


Potential  function  and  force  for  Cu  atoms  interacting 
through  a  Horse  potential  Vfrjr  D  f'exp(-lo<('r-r0))~  l  wp 
with  0:  o.  5*12.  eV,  <=(  =  /  •  i  S'  ^  5  /4'1  ,  z  L.S'il  5  A  . 

The  first  nearest  neighbor  is  located  at  the  minimum 
of  the  potential  (V^).  The  cut-off  range  is  1.6  t~ 
and  further  interaction  was  neglected. 
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(4) 

The  parameters  D,  o(  and  rQ  were  determined  byCotterill  and 
Girifalco using  experimental  values  for  the  energy  of  vaporization, 
the  compressibility  and  the  lattice  constant  of  f.c.c.  materials. 

Fixed  boundary  was  used  in  all  the  simulations.  Under  these 
conditions,  particles  on  the  boundaries  are  displaced  according  to  a 
continuum  elasticity  solution  and  remain  fixed  during  the  rest  of  the 
computation  (see  Chapter  Four). 

The  computer  simulations  were  carried  out  with  a  Morse  potential 
corresponding  to  Cu^  and  two  other  potentials  of  the  Morse  type  with 
the  parameter  o(  modified  (see  Table  5-1  at  the  end  of  this  chapter) . 
The  purpose  of  choosing  different  4  values  was  to  study  the  influence 
of  this  parameter  on  the  emission  of  dislocation  from  the  crack  tip. 

As  we  will  see  at  the  end  of  the  chapter,  this  parameter  seems  to 
affect  also  the  interaction  of  the  crack  tip  with  the  boundary,  the 
interaction  being  smaller  for  smaller  c<  . 

In  addition  to  the  CMD  experiments  to  determine  the  elastic  con¬ 
stants  of  the  system  (Chapter  4),  two  different  experiments  were  done 
to  find  the  ideal  cohesive  stress  and  the  ideal  shear  stress  <T-%  . 
The  ideal  cohesive  stress,  TCc  ,  can  be  determined  by  applying  a 
homogeneous  expansion  to  a  perfect  lattice  until  the  equally  extended 
bonds  reach  the  maximum  permissible  length  given  by  the  cut-off  of  the 
potential.  The  maximum  stress  reached  during  the  experiment  corres¬ 
ponds  to  the  ideal  cohesive  stress  CT^C  (Fig.  5.3).  The  ideal  shear 
stress,  Cn^  ,  can  be  determined  following  the  classical  experiment  of 
Frenkel ^  which  consists  of  the  shearing  of  two  rows  of  atoms  in  a 
homogeneously  strained  crystal  (Fig.  5.4).  As  noted  by  Orowan  the 
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Figure  5.3  Biaxial  stress  versus  strain  for  a  uniform  expansion 

of  a  perfect  lattice  whose  atoms  interact  through  a 
Morse  potential  (Cu) .  Point  a  on  the  curve  corresponds 
to  the  ideal  cohesive  stress,  .  Point  b  is  the 
rupture  point,  where  all  the  bonds  reach  the  cut-off 
distance  of  the  potential  function. 


DisPLflcerENT^ftTomc  spacepent 

Shear  stress  versus  displacement  for  the  shearing 
of  two  rows  of  atoms  which  interact  through  a  Morse 
potential  (Cu) .  Curve  a  was  obtained  with  displace 
ments  along  the  x  direction.  Curve  b  permitted 
simultaneous  relaxation  along  the  y  direction. 
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ideal  shear  stress  can  be  reduced  (curve  b)  if  atom  A  slides  over  atom 

C  instead  of  suffering  a  straight  displacement  (curve  a).  This  is  due 

to  the  steeply  repulsive  part  of  interatomic  potentials  generated  by 

the  overlapping  of  closed  shells. 

The  Burgers  vector  corresponding  to  a  certain  sheared  region 

.  (3) 

was  determined  by  the  equation 

l  =  l  JLl  (5*6) 

b  l 


where  the  derivative  was  approximated  using  the  differences  between 
the  displacements  of  two  nearest  neighbors. 

The  strain  fields  and  rotations  around  the  crack  tip  were  de- 

(Q\ 

termined  using  the  following  equation^  (Fig.  5.5): 


■xx 


'bu.  - 

a* 


■  .  ^  v  .  Ya  -  v, 

■  — 


4.  >V  .  U.Q-  v&  -^0  (  5  -  7  ) 

by  Ay  ax 

^xy  :  ^  -  c)  v  -  u.  a  -  U.0  ^  V a  -  Vo 

by  bx  Ay  Ax 


As  shown  in  Fig.  5.5,  the  relative  displacements  u  and  v  were 
interpolated  for  points  A  and  B  on  the  intersection  of  lines  between 
centers  of  the  nearest  neighbors  particles  and  the  positive  x  and  y 
directions.  Using  these  relative  displacements  and  the  initial 
separations  AXan^  Ay  ,  local  shear  strains,  dilatations 

and  rotations  were  determined.  With  the  strains, 
and  f  the  two  principal  strains  were  calculated.  These  results, 

*7 


97 


Y 


Fig.  5. 5  Interpolation  scheme  of  relative  displacements  of  bubbles  used 
in  computing  local  shears  and  dilatations  observed  in  sheared 
bubble  rafts. 
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jointly  with  the  rotations,  u/„v  ,  are  shown  in  Figs.  5.6  to  5.11  for 

7 

different  configurations  (before  and  after  relaxation)  and  three 
potential  functions.  Fig.  5.13. 

Figures  5.6a  and  5.6b  show  the  initial  and  final  configurations 
in  the  case  of  a  Morse  potential  with  «  equal  to  0.25,  M-3  potential. 

KTC  criterion  predicts  ductile  behavior  and  RT  criterion  brittle  be¬ 
havior  (Table  5.2).  Better  estimates  using  RT  criterion  can  be  made 

by  assigning  ,  the  dislocation  core,  more  realistic  values  deter- 
(3) 

mined  by  CMD. 

Figure  5.6b  shows  the  Burgers  vector  corresponding  to  the  sheared 
region  enclosed  in  the  Burgers  circuit  of  Fig.  5.12.  The  Burgers 
circuit  was  taken  along  a  region  of  perfect  material  and  finished  at 
point  b  to  add  up  the  effects  of  the  small  Burgers  vectors  correspond¬ 
ing  to  shear  along  the  planes  between  points  a  and  b.  The  total 
Burgers  vector  includes  also  the  effects  of  the  crack  opening  when  the 
linear  elasticity  theory  is  applied  to  generate  the  initial 
configuration. 

The  sheared  region,  when  a  large  stress  intensity  factor  is 
applied  (Fig.  5.6b),  extends  until  it  interacts  with  the  fixed  boundary. 
Figures  5.9a  and  5.9b  show  the  initial  and  final  configuration  for  the 
same  material,  M-3,  and  a  lower  stress  intensity  factor.  In  this  case 
the  sheared  region  does  not  interact  with  the  boundaries. 

Figure  5.6b  suggests  the  possibility  of  a  phase  transition  from 
a  triangular  to  a  cubic  system  by  rotation  of  the  lattice.  To  study 
this  kind  of  deformation  the  strain  field  and  rotations  were  deter¬ 
mined  using  the  Eq.  (5-7)  (Fig.  5.7a-b  and  Fig.  5.8a-b).  The  same 
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studies  were  done  with  a  different  Morse  potential,  M-2,  and  smaller 
shear  deformations  were  observed  (Fig.  5.10a-b).  The  different  be¬ 
havior  in  the  case  of  M-2  potential  seems  to  be  due  to  the  steeper 
slope  of  the  potential  function  (*<=0.4).  In  this  case,  the  inter¬ 
action  between  particles  has  changed  from  "soft  particles"  (HO. 25) 
to  "hard  particles"  (a  *0.4)  and  slip  along  certain  planes  is  now 
hindered  by  the  rigid  boundary  (i.e.,  dd'  Fig.  5.10b). 

By  increasing  the  parameter  ot  to  1.3,  M-l  potential, which  cor¬ 
responds  to  Cu  atoms,  the  interaction  with  the  boundary  is  even  greater 
(see  Fig.  5.11)  and  presumably  one  has  to  use  much  larger  systems  in 
order  to  observe  the  same  type  of  sheared  deformation. 

The  advantage  of  using  a  soft  potential,  M-3,  despite  the  un¬ 
realistic  value  of  o(  ,  is  that  the  weak  interaction  with  the  boundary 
may  enable  one  to  observe,  through  scaling,  the  behavior  of  a  more 
realistic  potential  such  as  M-l  in  a  larger  system. 

A  graphic  analogy  of  this  case  is  shown  in  the  bubble-raft 
(4) 

model  where  the  dislocation  width  changes  from  a  few  interatomic 
spacings  for  soft  bubbles  (1.9  mm  diameter)  to  approximately  40  inter¬ 
atomic  spacings  for  hard  bubbles  (0.3  mm  diameter).  This  explains 
again  the  difficulty  to  accommodate  the  dislocations  emitted  by  the 
crack  tip  in  a  few  interatomic  spacings  for  the  case  of  hard  potentials. 

KTC  criterion  predicts  ductile  behavior  for  M-2  and  M-3  poten¬ 
tials  .  However,  we  have  seen 

that  for  the  M-2  potential  the  interaction  with  the  fixed  boundary 
starts  to  affect  seriously  the  slip  along  some  planes.  In  the  case 
of  the  M-l  potential,  /  0T-t 


is  0.38,  too  close  to  0.5  to  make  a 
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clear-cut  prediction.  The  CMD  simulation  has  shown  crack  propagation 
by  breaking  bonds  along  the  crack  front  but  the  results  are  not  com¬ 
pletely  free  from  boundary  dependent  effects  by  the  reasons  given 
above. 

RT  criterion  seems  not  to  be  successful  in  predicting  the  cor¬ 
rect  behavior  for  potentials  M-2  and  M-3,  where  KTC  criterion  and  CMD 
both  predict  ductile  behavior.  As  noted  by  Rice  and  Thomson  the  ap¬ 
plication  of  linear  continuum  elasticity  in  the  region  close  to  the 
crack  tip  is  an  approximation  and  the  core  of  the  dislocation,  0  , 
has  been  estimated  according  to  the  Peiersls  model. 

A  comparative  study  of  the  prediction  of  these  two  criteria 
(9) 

has  been  made  by  Tyson  by  compiling  the  results  obtained  with  CMD. 
This  study  also  shows  that  KTC  criterion  is  more  successful 
than  the  RT  approximation. 

Finally,  the  Lennard-Jones  system  studied  in  Chapter  Four  clearly 
behaves  as  a  brittle  material.  Figs.  4.5a  to  4.5c.  This  behavior  is 
consistent  with  the  high  values  found  for  (Jc,  /  (T3t  and  §  ■ 
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Figure  5.6a  Initial  atomic  configuration  based  on  linear 

elasticity  theory  of  a  two-dimensional  crack 
embedded  in  an  infinite  medium.  System  consists 
of  436  particles  arranged  in  a  triangular 
lattice,  interacting  through  a  Morse 
potential  (M-3). 


5.6b  Same  as  Fig.  5.6a  after  the  relaxation 

of  the  atomic  region  (fixed  boundary) . 
Potential  M-3. 
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Figure 
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,7a  Principal  strains  at  the  initial  configuration 

(linear  elasticity  theory)  f  dilatation  strain, 
I  compressive  strain. 
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Figure  5.7b  Principal  strains  in  the  final  relaxed  con¬ 

figuration  (fixed  boundary)  f  dilatation 
strain,  |  compressive  strain. 
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Figure  5.8a  Rotations  at  the  initial  configuration 

(linear  elasticity  theory),  t  indicates 
the  magnitude  of  the  rotation. 
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Figure  5.8b  Rotations  at  the  final  relaxed  configuration 

(fixed  boundary) .  f  indicates  the  magnitude 
of  the  rotation. 
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Figure  5.9a 


Initial  atomic  configuration  (linear 
elasticity  theory)  Potential  M-3.  K- 0.5 
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Final  relaxed  configuration  (fixed  boundary) . 
Particles  interact  through  a  Morse  potential 
(M-3) .  Ks  0.5  Kc 


Figure  5.9b 
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Initial  configuration  (linear  elasticity 
theory) .  Potential  M-2.  K  -  Kc 


Figure  5.10a 


Figure 


.10b  Final  relaxed  configuration  (fixed 

boundary).  Potential  M-2. 
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Atomic  configuration  during  the  application 
of  flexible  boundary.  The  critical  stress 
has  been  exceeded  and  the  crack  propagates 
breaking  bonds.  Potential  M-l. 


Figure  5.11 
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Figure  5.12  Burgers  circuit  in  a  perfect  reference 

crystal.  The  circuit  has  been  finished 
at  point  b  to  add  all  the  Burgers  vectors 
corresponding  to  shear  on  planes  between 
points  a  and  b.  The  same  circuit  has  been 
applied  to  the  final  configuration  (Fig.  5.6b) 
to  determine  the  total  Burgers  vector. 
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Characteristic  parameters  and  properties  of 
several  materials  studied  in  CMD  simulation. 
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Table  5.2  Predicted  behavior  according  to  Kelly,  Tyson  and  Cottrell  (KTC) 

criterion.  Rice  and  Thomson  (RT)  criterion,  and  observed  behavior 
by  Computer  Molecular  Dynamics  simulation  (CMD). 
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Chapter  Six 


Summary  and  Conclusions 
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6.  Summary  and  Conclusions 

Computer  Molecular  Dynamics  has  been  applied  to  the  study  of 
fracture  in  a  two-dimensional  lattice.  The  major  purpose  of  this  work 
was  to  study  and  establish  interrelations  between  several  macroscopic 
properties  such  as  fracture  toughness,  stability  of  the  crack  against 
emission  of  dislocation  and  deformation  field  around  the  crack  tip. 

We  have  developed  a  computer  simulation  program  which  can  be  used 
with  crystalline  solids,  liquids  and  gases,  under  a  variety  of  boundary 
conditions,  periodic,  free,  fixed  and  flexible  boundaries. 

The  first  step  in  the  development  of  the  computer  simulation 
model  consisted  of  setting  up  the  standard  CMD  techniques  which  per¬ 
mitted  an  accurate  and  efficient  determination  of  the  classical  trajec¬ 
tories  of  the  particles  and  the  static  and  dynamic  properties  of  the 
system.  The  simulations  were  carried  out  in  a  perfect  monoatomic 
system  (rare  gas).  Additional  tests  were  done  to  study  the  accuracy 
of  the  calculations:  conservation  of  energy  and  reversibility  of 
trajectories . 

The  microscopic  elastic  constants,  dispersion  relation  and  phonon 
spectrum  of  the  system  were  determined  by  lattice  dynamics.  These 
calculations  were  not  necessary  for  the  rest  of  the  work  but  of  in¬ 
terest  in  order  to  study  anharmonic  effects  and  derive  analytical  ex¬ 
pressions  for  calculating  the  macroscopic  elastic  constants  of  the 
system. 

Two  boundary  conditions  have  been  studied  in  detail.  A  fixed 
boundary  condition  consists  of  placing  the  boundary  particles  according 
to  a  continuum  elasticity  solution.  A  flexible  boundary  condition 
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consists  of  surrounding  the  simulation  region  by  a  continuum  region 
in  which  particles  are  moved  according  to  a  Green's  function  to  relax 
the  forces  arising  from  the  interactions  between  the  atomic  system 
and  the  continuum.  The  Green's  function  used  is  that  derived  by  Hirth 
for  the  case  of  isotropic  elastic  media  with  a  crack,  a  result  we  have 
rederived  for  the  case  of  a  line  force  in  an  infinite  medium  following 
a  procedure  similar  to  the  determination  of  potentials  in  electro¬ 
statics.  The  simplicity  and  generality  of  this  method  permit  direct 
applications  to  the  determination  of  the  displacement  fields  of  other 
point  and  line  defects  (i.e.,  interstitials  and  dislocations). 

The  flexible  boundary  procedure  was  found  to  be  efficient  and 
accurate  in  studying  crack  tip  configurations  in  brittle  materials. 

It  was  possible  to  determine  the  critical  stress  values  for  propagating 
the  crack  with  a  smaller  number  of  particles  compared  to  the  fixed 
boundary  method. 

The  critical  stress  intensity  factor,  determined  by  CMD  simula¬ 
tion  in  brittle  materials,  has  been  found  to  be  greater  than  the  value 
predicted  by  Griffith's  theory.  However,  this  value  is  quite  sensitive 
to  the  assumed  crack  surface  interactions  in  the  simulation.  Higher 
values,  about  two  times  Griffith's  predictions,  were  found  when  inter¬ 
actions  were  permitted  between  atoms  situated  on  opposite  crack  sur¬ 
faces,  and  lower  values,  about  15%  greater  than  the  Griffith  criterion, 
were  found  when  no  interactions  were  assumed. 

The  difference  between  these  results  and  Griffith's  theory  can 
be  explained  as  follows.  First,  the  Griffith  criterion  is  based  on  an 
energy  balance  and  it  is  considered  a  necessary  but  not  sufficient 
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condition  for  crack  propagation.  Secondly,  the  equations  used  by 
Griffith  are  based  upon  macroscopic  linear  elastic  continuum  properties 
and  do  not  consider  the  local  atomic  relaxations  around  the  crack  tip 
and  the  changes  in  stress  distribution  with  different  crack  tip  atomic 
models.  In  the  case  of  plasticity  we  should  include  as  another  reason 
for  the  large  critical  stress  found  by  simulation  the  additional  energy 
expended  in  the  emission  of  dislocations.  In  the  test  of  Griffith’s 
criterion  we  have  worked  with  a  brittle  solid,  a  Lennard-Jones  rare 
gas  system.  Other  potentials  appropriate  for  brittle  materials  or  even 
truncated  Hookean  potentials  would  be  interesting  cases  for  study  by 
CMD. 

The  flexible  boundary  method  was  found  to  be  not  so  efficient, 
in  the  sense  of  being  applicable  to  systems  with  small  number  of  par¬ 
ticles,  when  high  stresses  were  applied  and  large  non-linear  regions 
were  produced  around  the  crack  tip.  The  constraint  imposed  by  the 
continuum  solution  at  the  boundaries  does  not  permit  one  to  study 
plastic  relaxation  by  the  emission  of  dislocation  from  the  crack  tip. 
However,  this  constraint  was  not  so  important  during  the  first  stages 
of  the  plastic  region  formation  or  when  small  stresses  were  applied. 

In  these  cases,  the  plastic  region  was  small  and  did  not  interact 
with  the  boundaries. 

The  strain  and  rotation  fields  were  determined  for  several 
stresses  and  different  potential  functions.  The  purpose  of  using  dif¬ 
ferent  potential  functions  was  to  study  the  influence  of  the  poten¬ 
tial  parameters  on  the  brittle  or  ductile  behavior  observed  by  com¬ 


puter  simulation. 
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In  the  case  of  ductile  materials,  which  correspond  to  Morse 
potentials  M-2  and  M-3,  emission  of  dislocation  from  the  crack  tip  was 
observed  and  the  size  of  the  sheared  region  increased  with  higher 
applied  stresses.  In  some  of  the  results  obtained  with  the  potential 
M-3,  analysis  of  the  rotation  field  around  the  sheared  region  suggests 
a  phase  transition  from  triangular  to  cubic  lattice. 

The  Kelly,  Tyson  and  Cottrell  (KTC)  and  Rice- Thomson  (RT) 
criteria  for  predicting  brittle  or  ductile  behavior  in  crystalline 
materials  were  studied  and  their  predictions  confronted  to  CMD  obser¬ 
vation.  Both  criteria  are  based  on  continuum  linear  elastic  deriva¬ 
tions  and  in  order  to  determine  their  predictions  the  elastic  con¬ 
stants,  surface  energy  and  ideal  stresses  were  calculated  by  indepen¬ 
dent  simulations  for  all  the  potential  functions. 

The  predictions  of  KTC  criterion  gave  good  agreement  with  CMD 
observations.  In  the  case  of  the  Morse  potential  M-l,  which  is  more 
steeply  repulsive  compared  to  M-2  and  M-3,  the  behavior  is  not  well- 
defined.  While  the  value  of  the  ratio  0c*  /  tr7c  is  close  to  0.5,  the 
limit  for  ductile  behavior,  brittle  behavior  was  observed  in  CMD  re¬ 
sults.  However,  the  computer  results  also  indicated  strong  inter¬ 
action  of  shear  deformation  with  the  boundaries  and  this  could  have 
been  the  reason  for  the  absence  of  dislocations  in  the  data. 

The  RT  criterion  seems  to  underestimate  the  possibility  of 
ductile  behavior.  However,  more  accurate  predictions  can  be  obtained 
by  using  CMD  to  determine  more  realistic  values  for  the  core  of  the 
dislocation.  The  core  width,  as  it  has  already  been  discussed,  can 
vary  significantly  with  the  slope  of  the  potential  function,  the 
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width  increasing  with  a  steeper  slope.  This  fact  would  explain  again 
the  absence  of  dislocations  in  the  simulation  with  the  M-l  potential, 
where  the  emission  of  dislocations  would  be  prevented  by  the  strong 
interaction  of  the  wide  and  uncomp ressible  core  of  the  dislocation  with 
the  boundaries. 

The  brittleness  of  the  material  can  be  affected  by  the  depth  of 
the  potential  function,  assuming  no  significant  changes  occur  to  other 
physical  constants.  The  surface  energy,  or  energy  necessary  to  create 
a  new  surface  by  breaking  bonds  is  logically  connected  with  this 
parameter  that  determines  the  energy  to  break  one  atomic  bond. 

Further  improvements  in  the  study  of  emission  and  propagation 
of  dislocations  from  the  crack  tip  could  be  done  by  determining  the 
Peirls  stress  and  core  width  of  the  dislocation  by  CMD,  using  the 
potential  functions  studied  in  this  thesis.  The  distance  travelled 
by  a  dislocation  in  the  presence  of  the  crack  tip  stress  field  will 
depend  ultimately  on  the  stress  necessary  to  overcome  the  friction  of 
the  lattice  (Peierls  stress),  and  it  will  define  the  extension  of  the 
plastic  region.  The  width  of  the  dislocation  core  would  help  to  pre¬ 
dict  accurately  the  behavior  of  the  system  using  the  RT  criterion  and 
also  would  explain  the  influence  of  the  boundaries  on  the  emission  of 
dislocations  from  the  crack  tip.  Flexible  boundary  conditions  could 
play  an  important  role  in  the  determination  of  the  Peierls  stress  and 
core  width  of  the  dislocation.  In  this  case  our  system  would  consist 
of  a  dislocation  core,  with  a  width  of  a  few  interatomic  spacings,  and 
an  exterior  linear  elastic  region  that  can  be  accurately  treated  with 
our  flexible  boundary  procedure. 
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Another  interesting  area  of  future  work  is  the  dynamics  of 
crack  propagation.  This  problem  is  well  suited  to  studies  by  CM), 
but  involves  difficulties  with  regard  to  appropriate  boundary 
conditions. 

Finally,  it  is  necessary  to  remark  that  all  the  present  simula¬ 
tions  were  carried  out  with  a  two-dimensional  system  at  0°K.  There 
are  known  physical  properties  such  as  the  ductile  behavior,  of  Cu 
which  may  appear  to  be  inconsistent  with  our  simulation,  in  this  case 
the  results  on  Morse  potential  M-l.  We  do  not  think  this  is  a  real 
conflict  because  there  are  close  packet  planes  in  a  three-dimensional 
system  which  are  not  present  in  a  two-dimensional  system.  It  would  be 
desirable  to  verify  directly  by  simulation  that  a  three-dimensional 
solid  with  potential  M-l  is  indeed  ductile.  The  consideration  of  this 
factor  would  not  involve  great  computational  difficulties  and  could 


be  another  area  of  future  work. 
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